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ABSTRACT 

We study the growth of the explosion energy after shock revival in neutrino-driven explosions 
in two and three dimensions (2D/3D) using multi-group neutrino hydrodynamics simulations 
of an 11.2M© star. The 3D model shows a faster and steadier growth of the explosion energy 
and already shows signs of subsiding accretion after one second. By contrast, the growth of 
the explosion energy in 2D is unsteady, and accretion lasts for several seconds as confirmed 
by additional long-time simulations of stars of similar masses. Appreciable explosion energies 
can still be reached, albeit at the expense of rather high neutron star masses. In 2D, the bind¬ 
ing energy at the gain radius is larger because the strong excitation of downward-propagating 
g-modes removes energy from the freshly accreted material in the downflows. Consequently, 
the mass outflow rate is considerably lower in 2D than in 3D. This is only partially com¬ 
pensated by additional heating by outward-propagating acoustic waves in 2D. Moreover, the 
mass outflow rate in 2D is reduced because much of the neutrino energy deposition occurs in 
downflows or bubbles confined by secondary shocks without driving outflows. Episodic con¬ 
striction of outflows and vertical mixing of colder shocked material and hot, neutrino-heated 
ejecta due to Rayleigh-Taylor instability further hamper the growth of the explosion energy in 
2D. Further simulations will be necessary to determine whether these effects are generic over 
a wider range of supernova progenitors. 

Key words: supernovae: general - hydrodynamics - instabilities - neutrinos - radiative trans¬ 
fer 


1 INTRODUCTION 


After decades of research, the mechanism powering core-collapse 
supernovae is still one of the outstanding problems in theoreti¬ 
cal astrophysics. The so-called delayed neutrino-driven mechanism 
currently remains the most popular and best explored explanation 
(see [Janka 2012, Burrows 2013 for current reviews) for “ordi¬ 
nary” supernova explosions not exceeding ~10 51 erg. In its mod¬ 
em guise, the neutrino-driven mechanism relies on additional sup¬ 
port for shock expansion in the form of hydrodynamic instabilities 


1992 Herant et al. 1994 Burrows, Hayes & Fryxell 1995|[Janka 

& Miiller[l996, Muller & Janka 1997) and the standing accretion 

shock instability (SASI, 

Blondin, Mezzacappa & DeMarino 2003, 

Laming 2007 Foglizzo 

et al. 2007, Guilet & Foglizzo 2012). In- 

deed, many successful n 
driven shock revival (m 
axisymmetry) have strei 
driven mechanism over t 

lulti-dimensional simulations of neutrino- 
ostly in 2D, i.e. under the assumption of 
ngthened our confidence in the neutrino- 

he recent years (Burns et al. 

2006a, Marek 

& Janka 20091 Yakunin et al. 2010 Suwa et al. 2010; Muller, Janka 

& Marek 2012, Muller, Janka & Heger 2012, Janka et al. 2012, 
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|Bruenn et al.|2013||Suwa et al.|2013[|Pan et al.|2015| >. However, 
both the long-time evolution of the successful 2D models as well as 
the advent of three-dimensional (3D) simulations have revealed two 
serious challenges for the neutrino-driven paradigm: With the ex¬ 
ception of |Bruenn et al.|(20T3) , the majority of successful 2D sim¬ 
ulations tended to produce explosions that are probably underener- 
getic. Moreover, the most ambitious self-consistent 3D simulations 
with multi-group neutrino transport have so far either failed to yield 
explosions ( [Hanke et al.|2013| |Tamborra et al.|[2014a|b| >, or, in the 
few successful cases, showed a considerable delay in shock revival 
(Mel son et al.|2015[ Lentz et al.|2015|) and signific antly smaller ex- 
plosion energies jTakiwaki, Kotake & Suwa|2014|>. Only the explo- 
sion of a 9.6M 0 star recently simulated by [Melson, Janka & Marek| 
( |2015| > is an exception from this trend. Whether the neutrino-driven 
mechanism provides a robust explanation for shock revival and can 
account for the observed explosion properties of core-collapse su¬ 
pernovae may appear doubtful in the light of these results. 


There is now an emerging consensus about the reasons un¬ 
derlying the more fundamental problem of missing or delayed ex¬ 
plosions in 3D. Both simple light-bulb and leakage-based simula¬ 
tions { [Hanke et al.|2012[|Couch|2013b|a[|Couch & Q’Connor|2014[ 
|Couch & Qtt|2015) as well as models with multi-group transport 
|Hanke et al.||2013[ |Takiwaki, Kotake & Suwaj2014| find an ar- 
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Figure 1. Density p and entropy s for the four progenitors si 1.0 (red), si 1.2 
(black), si 1.4 (light brown), and si 1.6 (green) as a function of enclosed 
mass. 


tificial accumulation of turbulent kinetic energy on large spatial 
scales in 2D due to the action of the inverse turbulent cascade 
( |Kraichnan| 1967) . Effectively, the forward cascade in 3D provides 
for more efficient damping/dissipation of the turbulent motions in 
the post-shock region, resulting in smaller turbulent kinetic ener¬ 
gies. Since the turbulent kinetic energy is directly related to the 
Reynolds stresses (i.e. loosely speaking, the “turbulent pressure”) 
that have been identified as the primary agent fostering shock ex¬ 
pansion in multi-D ( Burrows, Hayes & Fryxell|199 5 Murph y, Do-| 


|lence & Burro ws|2013[[Couch & Qtt|2015|lMiiller & Janka|2015[ ), 


this affects the critical neutrino luminosity L crit ( Burrows & Goshy| 

|1993[ |Murphy & Bu rrows 2008) required to power an explosive 
runaway. However, even though the higher explosion threshold in 
3D has emerged as a systematic effect, it nonetheless remains a 
small effect: Current light-bulb models (Hanke et al. 2012) |Couch| 
2013a, Dolence et al. 2013) invariably show a similar reduction of 
20... 30% in critical luminosity in 2D/3D compared to ID, with 
|Dolence et al.| ( |2013] ) still finding a slightly lower explosion thresh¬ 
old in 3D. Likewise, neutrino hydrodynamics simulation (Hanke 
|et al. ||2013l |Takiwaki, Kotake & Suwa|[2014| > show very similar 
heating conditions in 2D and 3D prior to shock revival (and even 
transient phases with better heating conditions in 3D in |Hanke et al.| 
|2013| ), although the small differences between 2D and 3D eventu¬ 
ally thwart shock revival in the 3D models of |Hanke et al.| ([2013 ) 
and Tamborra et al. (2014a b). Even though the adverse effects of 
the forward cascade in 3D may still be underestimated by the rela¬ 
tively crude grid resolution that current models can afford ( [Hanke| 
|et al.|2012||Couch|2013a[|Abdikamalov et al.|2014| ), it thus emerges 
that 3D models of neutrino-driven supernova explosions are very 
close to the explosion threshold. Consequently, relatively small re¬ 
finements in the simulations and the initial models may be suffi¬ 
cient to obtain explosions, e.g. moderate rotation ( [Nakamura et al.| 
|2014| ), magnetic fields ([Obergaulinger, Janka & Aloy 2014), as- 
phericities in the progenitor core ( [Couch & Ott|2013||Couch et al.| 
[2015, Muller & Janka 2015), or modifications to the standard neu¬ 
trino interaction rates ( [Melson et al. [20151 . The emergence of the 
spiral mode of the SASI could even allow strongly SASI-dominated 
models to explode earlier in 3D than in 2D ( |Fernandez|2015| >. 

By contrast, the problem of underenergetic neutrino-driven 
explosions in 2D has so far gone without a convincing theoret¬ 
ical explanation, and fewer suggestions have been made to rem¬ 
edy it, although it may be more serious in the sense that it affects 


even models with successful shock revival: While recent analy- 
of well-observed supernovae ([Tanaka et al.|T009| |Utrobin & 


Chug ai|2011| |Poznanski|2013[ Chugai & Utrobin|2014| Pejcha & 


Prieto 2015) suggest a broader range of explosion energies for type 
II-P supernovae between ~(0.1.. .2) x 10 51 erg instead of a sin¬ 
gle “canonical” value of 10 51 erg, there is arguably still a discrep¬ 
ancy, since almost none of the successful 2D and 3D models with 
multi-group neutrino transport show explosion energies consider¬ 
ably above 10 50 erg, e.g. t he final values at the en d of the simula¬ 
tions are a few 10 49 erg in Marek & Janka (2009) for proge nitors 
with 11.2M© and 15 M©, < 10 50 erg for a 13 M© progenitor in 


et al. 


Suwa 


< [2010] ), and 4 x 10 49 erg (11.2M© progenitor), 1.3 x 10 5U erg 


Janka et al. 


( 2012| ). More- 


(15M©), and 1.3 x 10 5 ° erg (27Af 0 ) in 
over, these estimates are not corrected for the “overburden”, i.e. the 
binding energy of the layers outside the shock, so that it remains 
unclear whether the explosions become energetic enough to shed 
the envelope at all. At first glance, the low explosion energies may 
simply be due to the fact that the simulations typically terminate 
before the explosion energy reaches its asymptotic value. While it 
can be argued that the final explosion energies may yet be higher 
by a factor of several because continued accretion sustains strong 
neutrino emission after shock revival that can power outflows from 
the gain region for >0.5 s, this assumption creates several other 
problems: Sustained accretion over >0.5 s might shift the result¬ 
ing remnant mass distribution well above the average birth mass of 
neutron stars M grav « 1.35M© ( |Schwab, Podsiadlowski & Rappa-| 
|port|2010| inferred from observations (which may, however, suffer 
from a selection bias). Only the 2D models of Bru enn et al.| ( |2014] ) 
form an exception from this trend; they obtain explosion energies 
in the range of (3.4... 8.8) x 10 5 ° erg for progenitors between 12M© 
and 25M© as well as reasonable Nickel masses. 


While this is encouraging, explosion energies above 10 51 erg 
still remain unexplained, and the problem of underenergetic super¬ 
nova explosion will arguably still linger as long as the discrepan¬ 
cies between the different simulation codes are not resolved. In¬ 
terestingly, [Melson, Janka & Marek| |2015| ) recently reported that 
3D effects, while apparently detrimental for shock revival in more 
massive progenitors, actually boost the explosion energy in a 9.6M© 
progenitor by ~10% by reducing the infall velocity in the accretion 
downflow and hence the cooling below the gain layer. Here, we 
further investigate their intriguing premise that 3D effects, while 
hurtful for shock revival, can prove beneficial in other situations 
and contribute to solving the problem of low explosion energies. 


In this paper, we present a successful 3D multi-group neu¬ 
trino hydrodynamics simulation of an 11.2M© progenitor with the 
CoCoNuT-FMT code ( [Muller, Janka & Dimmelmeier|2010[|Muller[ 
|& Janka|2015| ) as further evidence that 3D turbulence plays a posi¬ 
tive role after the onset of the explosion. By comparing the dynam¬ 
ics and energetics of the 3D explosion model to several long-time 
simulations of 2D progenitors in the mass range between 11M© 
and 11.6M©, we demonstrate how 3D effects can lead to a faster, 
more robust growth of the explosion energy provided that shock 
revival can be achieved. So far, the long-time effects of the dimen¬ 
sionality during the first hundreds of milliseconds to seconds after 
shock revival have received less attention than the role of the third 
dimension in shock revival: Successful first-principle models are 
still scarce and cannot be evolved for a sufficiently long time yet 
to address this phase in detail, while light-bulb-based studies of su¬ 
pernova energetics in 2D and 3D ( [Handy, Plewa & Odrzywofek| 
|2014[ ) cannot adequately account for the feedback of the subsiding 
accretion onto the neutrino heating and do not show the drawn-out 
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time [s] 


Figure 2. Shock propagation and diagnostic explosion energy £ e xpi for the 
11.2M© progenitor in 2D and 3D: The top panel shows the maximum, 
minimum (solid), and average (dashed) shock radius for model sll.2_2Da 
(black), s 11.2_2Db (blue), and s 11.2_3D (red). The middle and bottom panel 
show the diagnostic explosion energy £ e xpi and its time derivative dE cxp \/dt. 


long-time accretion characteristic of first-principle models. In this 
paper, we show that this phase deserves a closer look. 

Our paper is structured as follows: In Section [2] we review 
the numerical methods for hydrodynamics and neutrino transport 
used in our version of the CoCoNuT code, including a brief sketch 
of the modifications used in the 3D version. In Section [3] we first 
give a descriptive overview of the differences in shock propagation 
and explosion properties between the 2D and 3D models. Since the 
question of shock revival in 3D is not the objective of our current 
study, we only provide a brief, cautionary assessment of shock re¬ 
vival in the 3D model against the backdrop of recent first-principle 
models in Section [4] In Section [5] we then analyze the physical 
effects underlying these differences by combining a separate analy¬ 
sis of the outflows and downflows in the spirit of|Melson, Janka &| 
(M arek| < |20 1 5] ) . Finally, we summarize our results and discuss their 
implications in Section[6] 


2 NUMERICAL METHODS AND MODEL SETUP 
2.1 Progenitor Models 

We simulate the collapse and post-bounce phase of the (non¬ 
rotating) 11.2M 0 solar-metallicity progenitor si 1.2 of [Woosley,| 
|Heger & Weaver| < |2002| ) in 2D and 3D. In order to gauge the ef¬ 
fect of stochastic model variations, we perform two different 2D 



Figure 3. Shock propagation and diagnostic explosion energy Zs e xpi for the 
different progenitors in 2D: The top and middle panel show the maximum 
and average shock radius, respectively. The bottom panel shows the diag¬ 
nostics explosion energy £ exp i as a function of time (solid lines). Dashed 
lines show the time evolution £ exp i - E ow , i.e. the diagnostic energy cor¬ 
rected for the binding energy (overburden) E oy of the material ahead of the 
shock. Red, black, blue, light brown, and green curves are used for models 
sll.0_2D, sll.2_2Da, sll.2_2Db, sll.4_2D, sll.6_2D. 

simulations (sll.2_2Da and sll.2_2DtQ) for this progenitor, and 
we also conduct simulations for three other solar metallcity pro¬ 
genitors of |Woosley, Heger & Weaver] ( |20021 ) with similar zero- 
age main sequence (ZAMS) mass and density structure (11M©, 
11.4M 0 ,11.6M©). We randomly perturb the radial velocity v r at the 
beginning of collapse using a perturbation amplitude 6v r /v r = 10 -5 . 

In Figure [I] we show density and entropy profiles for the four 
progenitors simulated in the different 2D and 3D runs. Despite 
small differences in the location of the interfaces between the dif¬ 
ferent shells, the models are very similar in terms of structure with a 
pronounced density jump between the silicon shell and the convec- 

1 Scattering on nuclei (including cr-particles) was switched off after bounce 
for model sll.2_2Db, which leads to minor changes in early shock propa¬ 
gation, but is inconsequential for the long-time evolution. Since the energy 
exchange due to recoil in neutrino-nucleon scattering was taken to be pro¬ 
portional to the total scattering opacity instead of the neutrino-nucleon scat¬ 
tering opacity only (see Equation A31 in |Muller & Janka|2015) in model 
sll.2_2Da and all other models, the runs with neutrino scattering on nu¬ 
clei overestimate pre-heating from heavy flavor neutrinos outside the shock 
during the early post-bounce phase (whereas nuclei actually receive negli¬ 
gible recoil in scattering reactions), which leads to a slight reduction of the 
heavy flavor neutrino luminosity and a slightly slower contraction of the 
proto-neutron star compared to sll.2_2Db. 
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tive shell above the active oxygen burning zone. As we shall see, 
the 2D models of the different progenitors are qualitatively very 
similar in terms of explosion dynamics and energetics and should 
thus illustrate the generic behaviour of supernova explosions origi¬ 
nating from progenitors in this mass range. 


2.2 Numerical Methods 

The simulations are performed with the general relativistic (GR) 


neutrino hydrodynamics code CoCoNuT (Dimmelmeier, Font & 
Miiller|2002l |Miiller, Janka & Dimmelmeier|2010[|Muller &"j anka 


2015) . Our version of CoCoNuT uses piecewise parabolic (PPM) 
reconstruction ( [Colella & Woodwa rd|1984) combined with a hy¬ 
brid HLLC/HLLE Riemann solver < Mignone & Bodo|2005) to ob¬ 
tain higher-order spatial accuracy. CoCoNuT employs spherical po¬ 
lar coordinates (r, 6 , ip), which leads to strong time-step constraints 
near the polar axis in 3D due to the converging grid geometry. 
We circumvent this problem using an adapted version of the mesh 
coarsening scheme of |Cerda-Duran| ( [2009) : While the equations 
of hydrodynamics are solved on a fine grid with constant spac¬ 
ing Sip in longitude everywhere, a filter is applied to the solution 
after each time step to remove short wavelength noise in the in¬ 
direction by projecting the conserved variables onto piecewise lin- 
ear/quadratic[^] functions in “supercells” that contain 2 n{ - 6) fine cells 
in the ^-direction. The projection algorithm is implemented conser¬ 
vatively, and the slopes for the filtered solution are obtained using 
the monotonized-central (MC) limiter of | van Leer| ( |1977j ). The su¬ 
percell size 2 n(6,) is chosen such that nsinO > 1/2 is maintained 
at any latitude. This ensures that the allowed CFL time step at 
high latitudes is at most shorter by a factor of ~2 compared to the 
equatorial region, and limits the filtering to a region of 30° around 
the pole, which corresponds to 13.3% of the total volume. Simi¬ 
lar techniques have long been used in numerical meteorology, cp. 
|Kageyama & Sato| ( |2004) and Chapter 18 in |Boyd| ( [2001) . Polar 
filtering allows us to maintain the same effective angular resolution 
of 1.4° in 2D and 3D with grids of N r xN e = 550 x 128 zones (2D) 
and NrXNeXNy = 550 X 128 x 256 zones (3D, fine grid) covering 
the innermost 10 5 km of the star, respectively. 

Like any other solution to avoid the coordinate singularity and 
the excessive time step constraint near the axis such as Cartesian 
coordinates, overset grids ( Kageyama & Sato||2004[ |Wongwatha^| 


|narat, Hammer & Mtiller| 2010[ |Melson, Janka & Marek||2015) 

• cubed-sphere grids (Ronchi, Iacono & Paolucci 1996, Koldoba| 


|et al.||2002l |Zink, Schnetter & Tiglio||2008| Fragile et al.||2009) , 

this polar filtering procedure has specific advantages and disadvan¬ 
tages: Unlike Cartesian codes, polar filtering allows us to maintain 
spherical symmetry in the initial conditions and explicit symmetry¬ 
breaking terms can be avoided. Different from cubed-sphere grids, 
the grid remains orthogonal; and global conservation laws are eas¬ 
ier to enforce than for overlapping overset grids. On the other hand, 
projecting the solution to piecewise linear function effectively in¬ 
troduces an anisotropy in the numerical viscosity and diffusivity 
(an unwelcome effect that is minimized by overset or cubed-sphere 
grids but also manifests itself for Cartesian grids that are prone to 
the development of m - 4 modes). 

The space-time metric is computed using the extended con¬ 
formal flatness condition (xCFC JCordero-Carrion et al.|2009) . Be- 


2 For example, we use linear functions for the Eulerian density D and the 
mass fractions A; so that the conserved partial masses DXi are represented 
by quadratic functions. 



0 0.2 0.4 0.6 0.8 


time after bounce [s] 

Figure 7. Selected mass shell trajectories (black) for model sll.2_3D com¬ 
puted from spherically averaged density profiles. The trajectories start with 
roughly equal spacing in log r shortly before bounce. The plot also shows 
the maximum, average, and minimum shock radius (red), the gain radius 
(light brown), and the radii corresponding to densities of 10 11 g cm -3 and 
10 12 g cm -3 (blue). 


cause the asphericities in the gravitational field are small for non¬ 
rotating core-collapse supernovae we use the monopole approxi¬ 
mation for the gravitational field, i.e. the lapse function a, the con¬ 
formal factor 0, and the radial component of the shift vector only 
depend on r, and the non-radial components f3 e and of the shift 
vector are set to zero. 

For the neutrinos, we use the fast multi-group transport (FMT) 
scheme of Muller & Janka ( 2015), which is based on a stationary 
two-stream solution of the relativistic transfer equation that is com¬ 
bined with an analytic variable Eddington factor closure at low op¬ 
tical depths. This schemes includes general relativistic effects un¬ 
der the assumption of a stationary metric, but neglects velocity- 
dependent effects like Doppler shift and aberration. The neutrino 
rates include emission, absorption, and elastic scattering by nu¬ 
clei and free nucleons (along the lines of |Bruenn|1985) as well as 
an effective one-particle rate for nucleon-nucleon bremsstrahlung 
and an approximate treatment of the energy exchange in neutrino- 
nucleon scattering reactions. Comparisons of the FMT scheme with 
the more sophisticated relativistic neutrino transport solver Ver¬ 
tex ( [Rampp & Janka|2002[ |Muller, Janka & Dimmelmeier|2010) 
showed excellent qualitative and good quantitative agreement. For 
more details, we refer the reader to Muller & Janka |2015) . 

In order to further alleviate the time-step constraint, the in¬ 
nermost part of the computational domain (where densities ex¬ 
ceed ~5 x 10 11 g cm -3 ) is calculated in spherical symmetry using 
a conservative implementation of mixing-length theory for proto¬ 
neutron star convection, a procedure that has been used in the con¬ 
text of supernova simulations before (e.g. | Wilson & Mayle|1988[ 
|Hiidepohl|2014) . The transition density is adjusted such that it lies 
inside the convectively stable cooling layer. 

In the high-density regime, we use the equation of state (EoS) 
of |Lattimer & S westy (1991) with a bulk incompressibility modu¬ 
lus of nuclear matter of K - 220 MeV. At low densities, we employ 
an EoS accounting for photons, electrons and positrons of arbitrary 
degeneracy, an ideal gas contribution from baryons (nucleons, pro¬ 
tons, Qf-particles and 14 other nuclear species), Nuclear reactions 
are treated using “flashing” as described in Rampp & Janka ( 2002). 
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Figure 4. Specific entropy for model sll.2_2Da (top row) and model sll.2_3D (in a slice almost perpendicular to the equatorial plane, bottom row) at 
post-bounce times of 80 ms, 140 ms, and 181 ms (left to right). Note that a different color scale for the entropy is used for each of these snapshots. 


Figure 5. Specific entropy for model sll.2_2Da (top row) and model sll.2_3D (in a slice almost perpendicular to the equatorial plane, bottom row) at 
post-bounce times of 241 ms, 471 ms, and 944 ms (left to right). Note that a different color scale for the entropy is used for each of these snapshots. 
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6 B. Miiller 



Figure 6. Volume rendering of the entropy in model si 1.2_3D at post-bounce times of 89 ms (top left), 134 ms (top right), 210 ms (bottom left), and 580 ms 
(bottom right). 


Table 1. Overview of Simulations. The extrapolation of the final remnant masses (last column) is discussed in Section |T3| 


Model 

Progenitor 

Dimensionality 

Post-Bounce 
Time Reached [s] 

Diagnostic Energy 
Reached [erg] 

Baryonic Neutron Star 
Mass Reached [M 0 ] 

Extrapolated Baryonic 
Remnant Mass [M 0 ] 

sll.0_2D 

sll.O 

2 

8.195 

1.3 xlO 50 

1.62 

1.62 

sll.2_3D 

si 1.2 

3 

0.944 

1.3 x 10 50 

1.33 

1.48 

sll.2_2Da 

sll.2 

2 

1.044 

5.0 x 10 49 

1.37 

— 

sll.2_2Db 

sll.2 

2 

6.003 

7.8 x 10 49 

1.47 

1.69 

sll.4_2D 

sll.4 

2 

6.129 

1.0 x 10 5 ° 

1.56 

1.63 

sll.6_2D 

si 1.6 

2 

11.453 

2.1 x 10 5 ° 

1.62 

1.63 
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Figure 8. Ratio v s h/v s h,MM between the angle-averaged shock 
velocity v s h = dr s h, aV g/df and the shock velocity v s h,MM = 
0.794(F exp i /M exp i) 1 /2 [M exp i /(p pre r 3 )]°- 19 predicted by the model of 
Matzner & McKee (1999). Note that r s h, aV g and its numerical derivatives 
need to be smoothed considerably to allow for a useful comparison. 
v s h/%,MM is only shown until 5 s after bounce because the automatic 
smoothing procedure becomes ineffective toward the end of simulations 
sll.0_2D, sll.2_2Db, and sll.4_2D and v s h/v s h,MM becomes highly 
oscillatory. 



Figure 9. Angle-averaged, density-weighted velocity profiles for model 
sll.2_3D at different post-bounce times. At the end of the simulation, the 
angle-averaged velocity is positive outside a mass coordinate of 1.35Af 0 , 
but the zero point is still moving outward in mass. Note that the angle- 
average extends over the post-shock and pre-shock region and cannot be 
used to infer the post-shock velocity directly. 

3 OVERVIEW OF SIMULATION RESULTS 

In all our simulations, runaway shock expansion sets in when the 
Si/SiO interface reaches then shock and the mass accretion rate 
drops rapidly. Figures [2] (all 2D/3D 11.2M© models) and [^(long¬ 
time evolution of all 2D models) provide an overview over the prop¬ 
agation of the shock and the growth of the explosion energies for 
the different models; they show the maximum, minimum (only Fig¬ 
ure [2), and angled-averaged shock radius, as well as the “diagnos¬ 
tic explosion energy” E exp i, which we define as the total net energy 
(i.e. gravitational+internal+kinetic energy) of all the material that 
is nominally unbound and is moving outward with positive radial 
velocity at a given time (cp. |Muller, Janka & Marek|2012||Bruenn| 


Figure 10. Baryonic neutron star masses (comprising all matter at densi¬ 
ties higher than 10 11 g cm -3 ) for the different 2D and 3D simulations as a 
function of time. 

|et al.|2014) . The nucleon rest masses are not included in the internal 
energy, i.e. nucleon recombination only contributes to the diagnos¬ 
tic energy once it actually takes places. Figure[2]also shows the time 
derivative of the diagnostic energy. Key results of the simulations, 
including the diagnostics energy and the baryonic remnant mass at 
the end of the simulations as well as estimates for the final remnant 
mass (see Section |33| below), are given in Table[I] 

3.1 Differences Between 2D and 3D During the First Second 

For the 11.2M© progenitor, the first second after bounce is shown 
in detail in Figure [2] both in 2D and 3D. In addition, Figures [4] 
and [5] illustrate the multi-dimensional flow morphology for mod¬ 
els sll.2_2Da and sll.2_3D on meridional slices, and 3D ray-cast 
images of neutrino-heated convective bubbles in model sll.2_3D 
before and after shock revival are shown in Figure [6] 

Prior to the infall of the Si/SiO interface, we find very similar 
shock trajectories independent of dimensionality. However, prompt 
convection develops slightly differently in 2D and 3D, and its resid¬ 
ual effect on the entropy and lepton number profiles leads to a slight 
divergence between the 2D and 3D models already at early times 
in many quantities (neutron star radius, gain radius, cooling pro¬ 
files, etc.). This effect is not unphysical per se , but is most probably 
exaggerated in our models because the FMT scheme tends to over¬ 
estimate the strength of prompt convection. In view of the large 
systematic effects that we shall discuss later, it is also inconsequen¬ 
tial, but needs to be borne in mind when comparing the different 
models. 

After the infall of the Si/SiO interface, the shock expands 
slightly faster in 3D than in 2D, and the explosion energy starts to 
reach appreciable positive values several tens of milliseconds ear¬ 
lier. Snapshots of the entropy for models sll.2_3D and sll.2_2Da 
during this phase can be seen in the middle and right columns of 
Figure [4] which show the development of large convective plumes 
in both cases. The reader will note that the plumes are somewhat 
aligned with the coordinate axis in 3D, which is clearly a result 
of the coordinate choice but need not be considered harmful as dis¬ 
cussed in Section]?] At later times, the morphology of the 3D model 
is quite different; instead of the broad, laminar downflows charac¬ 
teristic for 2D explosions, the interface between the downflows and 
the hot, neutrino-heated ejecta eventually becomes turbulent dur- 
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ing the infall, resulting in corrugated downflows and partial mixing 
with the neutrino-heated ejecta, as can be seen most perspicuously 
in the middle column of Figure [5] 

Soon after shock revival, the 2D models start to go through 
episodes of halting shock expansion or even transient shock reces¬ 
sion. While the growth rate dE oxp \/dt of the diagnostic explosion 
energy reaches values comparable to 3D for 100... 200 ms, the ex¬ 
plosion energy grows much less steadily in the long term and has 
reached only (4... 5) x 10 49 erg after 1 s. By contrast, the 3D model 
shows a steady growth of the explosion energy (1.3x 10 5 ° erg by the 
end of the simulation), and considerably faster shock expansion. As 
illustrated by the mass shell trajectories in Figure[7] the spherically 
averaged radial velocity behind the shock becomes positives about 
300 ms after bounce, and the mass shells shocked later than 500 ms 
after bounce appear to move outward steadily instead of eventually 
falling back onto the proto-neutron stars. 


3.2 Shock Propagation During the First Seconds 


Before we attempt to extrapolate the final remnant masses, it is use¬ 
ful to point out a simple analytic relation between the diagnostic 
energy and the shock velocity. During the later phases of the explo¬ 
sion when the explosion energy has saturated, simple analytic mod- 
e ls (|Sedov|1959||K ompanee ts|1960||Laumbach & Probstein|1969| 


Klimishin & Gn atyk|198T |Koo & McKee|1990||Matzner & McKee 


1999) provide a useful qualitative description of shock propagation 
in hydrodynamical simulations ([Woosley & Weaver 1995||Kifoni-| 

|dis et al.|2003||Wongwathanarat, Muller & Janka|2015} . These un¬ 
derlying models typically rely on the assumption of self-similarity 
and/or exponential or power-law approximations for the envelope, 
neglect the effect of gravity, do not account for continuous energy 
input into the ejecta, and have been derived under the assumption 
of spherical symmetry. During the first seconds covered in our sim¬ 
ulations, all these conditions are violated. It is remarkable that the 
approximate formula of |Matzner & McKee] ( [1999} , 


v S h,MM - 0.794 


^expl / 

m \ 


Ppre' 


( 1 ) 


nonetheless provides a reasonable estimate for the shock velocity 
v sh already a few hundreds of milliseconds after shock revival if it 
is evaluated using appropriate definitions for the explosion energy 
.Eexpi? the “ejecta” mass m, and the pre-shock density p pre : We find 
that Equation ([lj works well if v sh is taken to be the angle-averaged 
shock velocity, i.e. the time-derivative of the angle-averaged shock 
radius r sh avg , if the pre-shock density is evaluated at r sh avg , if £ expl 
is identified with the time-dependent diagnostics energy, and if the 
“ejecta” mass includes the entire mass enclosed by the shock from 
above and the gain radius from below (and thus cannot properly be 
termed “ejecta” mass as in the original work of |Matzner & McKee] 
1999). This is illustrated in Figure [8] which shows the ratio of the 
angle-averaged shock velocity v sh and the analytic estimate v sh ,MM 
of Matzner & McKee (1999). While there is considerable scatter, 
the numerical models fall in a band with 1 < v sh / v sh ,MM < 1-6 
most of the time, especially at time later than 1 s after bounce. Our 
models suggest v sh = 1.3v sh ,MM as a good analytic estimate for early 
shock propagation in core-collapse supernovae. 


3.3 Explosion Energies and Neutron Star Masses 

Although the explosion energy in model sll.2_3D has not yet 
reached its final value and there is still some accretion onto the 



Figure 11. Comparison of the post-shock velocity computed from Equa¬ 
tion {TJ (red curve) to the required shock velocity ft / (J3 - l)v esc for the sep¬ 
aration of outgoing and infalling mass shells for two different values of the 
compression ratio p (black curves) for model si 1.2_3D. 


proto-neutron star, there is no doubt that the incipient explosion 
will eventually expel the envelope. This is not only clear from the 
steady outward movement of the shocked mass shells at late times 
visible in Figure [7] and in the velocity profiles depicted in Figure [9} 
the diagnostic explosion energy is also significantly higher than the 
residual binding energ y of the pre-shock matter ( the “overburden” 
in the terminology of Bruenn et al.||2013 2014} of 5 x 10 49 erg. 
This is clearly different from models si 1.2_2Da and si 1.2_2Db and 
the similar 2D explosion models of the same progenitor discussed 
in |Buras et al.| ( [2006a} ; |Marek & Janka] ( |2009} ; |Miiller, Janka &| 
|M are k M20l2) . 


Nonetheless, the long-time simulations of the 2D models over 
several seconds show that even such supposedly tepid models even¬ 
tually develop steady shock propagation and reach sufficiently high 
diagnostic explosion energies to shed the envelope (Figure [1}. For 
the low-mass progenitors simulated here, the outgoing and infalling 
mass shells inevitably separate as the post-shock velocities behind 
the entire shock front become positive once it reaches the edge of 
the relatively small C/O core (1.7... 1.8M 0 ), which happens al¬ 
ready after a few seconds in these models. The acceleration of the 
shock at the steep density gradient between the C/O core and the 
He shell and the small binding energy of the He shell then result 
in a steady outward movement of the shocked matter. At that point, 
the overburden of the unshocked envelope becomes almost negligi¬ 
ble (e.g. 10 49 erg for sll.2_2Db, 5 x 10 48 erg for sll.6_2D), and we 
can determine relatively firm lower limits for the final explosion 
energy. 

Continuous accretion over several seconds provides for suf¬ 
ficient neutrino heating to power outflows and pump additional 
energy into the ejecta over this long time-scale, albeit at a rather 
modest rate. As a result, models that appear woefully underener- 
getic during the first second can still develop appreciable explosion 
energies, the best example being the 11.6M 0 model, where E e x P i 
grows from 3.5 x 10 49 erg at 1 s to 2.0 x 10 5 ° erg after 11s. The ex¬ 
plosion energies obtained after several seconds are comparable to 
the 3D case and compatible with supernova explosion energies at 
the lower end of the observed spectrum (see, e.g.,|Pejcha & Prieto| 
|20l5} . 


The fact that the diagnostic explosion energies increase more 
or less steadily over several seconds in the 2D models (except for 
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transient phases where the explosion geometry changes because a 
neutrino-driven outflow is shut off as discussed in Section [5. 2. 4| > 
has important implications for the usefulness of the diagnostic en¬ 
ergy as a predictor of the final explosion properties. |Perego et al.| 
( |2015| ) have recently pointed out on the basis of artificial ID explo¬ 
sions of |Perego et al.| ( [20T5] > that the diagnostic energy overshoots 
the final explosion energy and only approaches its asymptotic value 
very slowly on a time-scale of seconds, and suggest that a better es¬ 
timate for the final explosion energy can be obtained by subtracting 
the overburden E ow , i.e. the binding energy of the mass shells out¬ 
side the shock, from £ exp i- As illustrated by the comparison of E exp i 
and E cxp] - E ow in Figure |3| our 2D models show a somewhat dif¬ 
ferent behaviour; similar to the simulations of |Bruenn et al.|j20T4| , 
there is no overshooting of E oxp[ above its prospective final value 
for which it appears to furnish a lower bound rather than an upper 
bound. This is the result of a fundamentally different way to power 
the explosion in multi-D compared to ID: Once an explosion is 
triggered in ID, the outflow rate quickly drops and only a weak 
neutrino-driven wind can still pump energy into the ejecta over 
time-scales of seconds. The accumulation of shocked material with 
negative total energy therefore quickly dominates the total energy 
budget of the ejecta region and £ exp i decreases, while E exp] - E ow 
remains roughly constant by virtue of total energy conservation. 
The case for E txp] - E ow as a more compelling predictor of the final 
explosion energy is weaker in multi-D, however, where neutrino- 
driven outflows can continuously pump energy into the ejecta at a 
high rate, and a considerable part of the shocked material with neg¬ 
ative total energy is channeled onto the proto-neutron star instead 
of being swept along by the ejecta and reducing the diagnostic en¬ 
ergy. £ e x P i may still decrease somewhat on time-scales longer than 
5 ... 10 s as the shock propagates through the helium shell, and this 
introduces a residual uncertainty of up to ~ 15% in the final explo¬ 
sion energy, which we expect to lie in the range bracketed by £ exp i 
and £ exp i - E oy . It is also noteworthy that £ exp i - E ow does not ap¬ 
pear to be a good predictor for the final explosion energy at early 
times simply because its rise phase is much more drawn out than 
in artificial ID explosions and it only becomes positive ~1 s after 
bounce or later. 

While our simulations reach final explosion energies of the 
order of 10 50 erg, this comes at the expense of rather high neutron 
star masses: Figure [lO] shows that the baryonic neutron star masses 
M by in the 2D models all end up at values > 1.47M© and will def¬ 
initely exceed 1.6M© in cases like sll.0_2D and sll.6_2D. Unless 
selection effects favor the production of less massive neutron stars 
in binary systems for some reason, this potentially presents a seri¬ 
ous conflict with the inferred neutron star mass distribution. Even 
if the lowest-mass neutron stars are presumed to originate from 
electron-capture supernovae, the masses of the neutron stars in the 
2D models would end up well above the mean value of inferred 
baryonic masses of 1.5M© ( |Schwab, Podsiadlowski & Rappaport] 
|2010| . Since the simulated models represent progenitors with rela¬ 
tively small cores and a relatively early onset of the explosion, this 
is highly problematic. It is interesting to note that the 2D models 
of |Bruenn et al.|(2014| l also show such a tendency towards high 
neutron star masses despite their relatively high explosion energies 
(although this tendency is less striking than in our long-time sim¬ 
ulations), with M by = 1.461M© for their 12M© model B12-WH07 
and values well above 1.6M© for the three remaining simulations. 

The faster rise of the explosion energy in 3D could help to re¬ 
solve this discrepancy. Although the neutron star mass has not yet 
converged to a final value, the spherically-averaged velocity pro¬ 
files (Figure[9]indicate that the final “mass cut” is slowly emerging. 


At the end of the simulation, the net mass accretion rate onto the 
neutron star in si 1.2_3D is lower by a factor of ~2 compared to the 
corresponding 2D models. 

To obtain a quantitative estimate for the final neutron star 
mass, we follow Marek & Janka ( 2009), who argued that accretion 
must subside once the post-shock velocity v post becomes compara¬ 
ble to the escape velocity v esc - For a strongly asymmetric explosion 
geometry, v post is of course strongly direction-dependent. Hence 
material ahead of the neutrino-heated plumes originating from a 
given mass coordinate m in the progenitor will be accelerated to a 
higher post-shock velocity by the shock than material with the same 
initial m that is hit later in a direction where the shock expands more 
slowly, so that the actual “mass cut” does not correspond to a single 
mass shell m in the progenitor. Instead, the dividing line in initial 
mass coordinate will depend on angle. Nonetheless, one can argue 
that the criterion v pos t = v esc still yields a fairly reliable estimate for 
the final mass of the neutron star if an appropriate spherical average 
for v P ost is used. 

Equation ([lj for the average shock velocity allows us to 
extrapolate the evolution of the v post if necessary to estimate a 
spherically-averaged “mass cuf’^Jlf the pre-shock velocity is as¬ 
sumed to be negligible, the post-shock velocity becomes 


v P ost — 


jg-1 

p 


V s h, 


( 2 ) 


in terms of the ratio p of the post- and pre-shock density, and equat¬ 
ing this to the escape velocity yields the criterion 


P — 1 /2G(M by + Mg abl ) 

■7“ Vsh= V- 7 -* < 3 > 

where we include the entire mass interior to the shock and not just 
the mass of the neutron star when computing the escape velocity. 
At late stages, the compression ratio p typically drops below the 
value p = 7 for a radiation-dominated ideal gas with adiabatic in¬ 
dex y — 4/3 because of nuclear burning and/or because the strong 
shock approximation is not strictly applicable over the downflows. 
We therefore compare v sh in model sll.2_3D to the critical veloc¬ 
ity PUP - l)v esc for two different values of p in Figure m Fig- 
ure mi suggests a final baryonic remnant mass of 1.41... 1.48M© 
for sll.2_3D (to which late-time fallback might be added). This 
would imply that the shock has already passed the initial mass cut 
in some directions. Estimates along the same lines for the long¬ 
time simulations in 2D yield baryonic remnant masses of 1.62M© 
for si 1.0_2D, 1.63M© for sll.4_2D, 1.69M© for sll.0_2Db, and 
1.63 M q for sll.6_2D assuming p - 4. 

Using the approximate formula of |Timmes, Woosley &| 
1996) for the gravitational neutron star mass M grav , 

M grav * M by - 0.075A4 , (4) 


Weaver 


which provides a reasonable fit across different nuclear equations of 
states, the estimated baryonic neutron star mass for s 11.2_3D can be 
converted to a gravitational mass of 1.34M©, which would be well 
within the range of observed neutron star masses and slightly below 
the mean value of the higher-mass population of neutron stars from 
iron core progenitors suggested by [Schwab, Podsiadlowski & Rap-| 
|paport| |2010| l. For the 2D simulations the estimated gravitational 


3 Equation {TJ is also more convenient to use from the numerical point of 
view because the computation of the shock velocity as a numerical deriva¬ 
tive of the shock position typically yields very noisy results. 
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time after bounce [s] 

Figure 12. Comparison of key properties of the neutrino-driven outflows in 
2D (black curves) and 3D (red curves) as measured at a radius of 400 km. 
The top panel shows the mass outflow rate M out . The middle panel shows 
the total enthalpy flux /\ ou t as defined in Equation J6| (thick lines) along¬ 
side the time derivative dE exv \/dt of the explosion energy (thin lines). The 
average total enthalpy h out (thick lines) and total energy e out in the outflows 
are shown in the bottom panel. 

masses are higher by more than 0.1 M©. The 3D effects responsi¬ 
ble for the steeper rise of the explosion energy thus improve the 
agreement with the observational constraints considerably. 


4 ASSESSMENT OF SHOCK REVIVAL IN THE 3D 
MODEL 

In the remaining part of the paper, our main thrust will be to ex¬ 
plain the physical mechanisms behind the pronounced differences 
between 2D and 3D simulations in the explosion phase presented in 
Section[3] We do not investigate the differences in the pre-explosion 
phase, because the numerical methodology used in this study only 
allows limited conclusions concerning the problem of shock revival 
in 3D for reasons detailed below. Nonetheless, a few remarks about 
shock revival in model sll.2_3D are in order, if only to motivate 
why the remainder of this paper focuses completely on the explo¬ 
sion phase, and why simulations with a more rigorous treatment of 
the neutrino transport and the neutrino rates are needed to decide 
the fate of this particular progenitor model. 

Superficially, our results for the 11.2M 0 star in 3D may ap¬ 
pear to be at odds with the recent core-collapse supernova sim¬ 
ulations with multi-group neutrino transport find either no shock 
revival at all in 3D or only delayed shock revival compared to 3D 
( |Hanke et al.|2012||Hanke|20T4l|Tamborra et al.|2014a||Lentz et al.| 


|2015||Melson et aL||2015|). Specific ally, the 11.2M© model failed 
to explode in 3D ( [Tamborra et al.||2014a| in a simu lation using 
the Vertex-Prometheus code ( |Rampp & Jank a 2002, Bur as et al.| 
|2006b) . However, one should not attach undue importance to the 
different outcomes of the Vertex-Prometheus and CoCoNuT-FMT 
models: Although reasonably close agreement with the more rig¬ 
orous transport scheme in Vertex can be reached with the FMT 
scheme, the differences noted by Muller & Janka ( 2015]) are suf¬ 
ficiently large to matter for a marginal model like si 1.2. The fact 
that we find an explosion merely underscores how close the 11.2M 0 
Vertex model of Hanke (2014) and Tamborra et al. ( 2014a) comes 
to an explosive runaway, and that the extreme sensitivity of the su¬ 
pernova problem to the neutrino transport treatment and the micro¬ 
physics ( [Lentz et ah 2012b a, Melson e t al.|2015] > requires highly 
accurate first-principle models in order to reliably decide whether 
an individual progenitor close to the explosion threshold explodes 
or fails (although “imperfect” simulations may already unearth 
much of the relevant physics from such cases). Considering that the 
comparison between the FMT scheme and Vertex revealed slightly 
better heating conditions for a 15M 0 progenitor at early times, and 
that the average shock radius initially expands somewhat further in 
3D than in 2D in the simulations of the 11.2M 0 progenitor with 
Vertex, the different outcome of the FMT and Vertex runs is by 
no means unexpected. 


Potentially, the emergence of large-scale bubbles aligned with 
the coordinate axis (Figures[4][5]and[6]) also helps in pushing the 3D 
model above the explosion threshold at an early time (cp. |Thomp-| 
|son|2000| |Dolence et al.|2013l|Fernandez|2015| for the role of the 
bubble size in the development of runaway shock expansion). The 
alignment is clearly a consequence of our coordinate choice and 
may also be connected to the coarsening procedure for the polar 
region. Such artifacts are unavoidable for standard spherical polar, 
cylindrical, or Cartesian coordinates (where they manifest them¬ 
selves as a preferred excitation of m = 4 modes instead), because 
the grid geometry and spacing dictates the effective numerical dif- 
fusivity and viscosity of a code, and physical instabilities will grow 
preferentially in directions where they are least suppressed (or even 
aided) by numerical dissipation. If there is sufficient time for insta¬ 
bilities like convection and the SASI to reach saturation and go 
through several overturn time-scales or oscillation periods, these 
initial artifacts from the growth phase are eventually washed out, 
but in a situation where the growth of certain modes accelerates 
rapidly (e.g. after the infall of the Si/SiO interface) and then freezes 
out they can subsist throughout the simulation. Nonetheless, we do 
not view this a concern; the convective flow does not show any grid 
alignment prior to the infall of the interface, and outflows eventu¬ 
ally develop in the equatorial plane as well. Moreover, we found 
no grid alignment of sloshing/spiral motions in SASI-dominated 
models (which will be reported elsewhere). In more realistic sim¬ 
ulations, the explosion geometry will be dictated by anisotropies 
in the initial model, e.g. due to convective nuclear burning (| Arnett 

T994l|Bazan & Arnett| 1 9941 |T998l |AsIda & Arnett|2000l |Kuhlen, 

Woosley & Glatzmaier|20 03, Meakin & Arnett 2006 2 007b|a||Ar^ 

nett & Meakin||2011[ |Couch & Qtt||2015| ) or rotation. In a sense, 

the alignment of the most prominent high-entropy bubbles with the 
axis in model sll.2_3D is even fortunate for our further analysis 
because it eliminates the unavoidable grid alignment of 2D explo¬ 
sion models as a potential cause for the different energetics in 2D 
and 3D. 
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Figure 13. Volume-integrated heating/cooling rates (Jheat and Qcoo\ in the 
gain and cooling region for models sll.2_2Da and sll.2_3D. The inner 
boundary of the cooling region is defined (somewhat arbitrarily) by a den¬ 
sity of 10 13 g cm -3 . Note that a different scale is used for both rates. 



time after bounce [s] 

Figure 14. Properties of the gain radius in 2D and 3D: The top panel shows 
the binding energy (i.e. the sum of the gravitational, kinetic, and internal 
energy) at the gain radius in 2D (black solid curve) and 3D (red solid curve) 
alongside the Newtonian potential of the neutron star GM/ r ga i n (dashed 
lines), which sets the typical energy scale at the gain radius. Here, M is the 
neutron star mass, and r ga i n is the gain radius. The estimate for e ga \ n from 
Equation Gl}. which is based on the assumption that the Bernoulli integral 
at the gain radius is zero, is shown in blue for comparison. The middle panel 
shows r ga i n itself, and the bottom panel shows the angle-averaged tempera¬ 
ture at the gain radius, T ga \ n . 



Figure 15. Outflow efficiency 77 out , for models sll.2_2Da and sll.2_3D. 
rj out is defined as the ratio between the actual mass outflow rate M out and a 
fiducial scale <2heat/kgainl for the mass loss rate, see Equation {I 4 } . Note 
that 77 ou t is not limited to values r] out ^ 1 because fresh matter for the 
neutrino-heated outflows is also supplied by lateral mixing with the down¬ 
flows above the gain radius (where the matter is less tightly bound than at 
the gain radius) and because recombination also partly contributes in lifting 
the material out of the gravitational potential well. 

5 ANALYSIS OF 2D/3D DIFFERENCES 

We now turn to the underlying physical mechanism responsible 
for the different evolution of the 2D and 3D models during the 
explosion phase. The first step towards understanding the differ¬ 
ent dynamics of the 2D and 3D models consists in considering the 
outflows and downflows separately (in the vein of |Melson, Janka| 
|& Ma rek 2015) to analyze the injection of mass and energy into 
the “ejecta region” with positive binding energy (similar to Bruenn 
|et al.||2014| >. We partition the computational domain into two re¬ 
gions with positive radial velocity (v r > 0) and negative radial ve¬ 
locity (v r < 0) and then compute total fluxes and averages of several 
hydrodynamic quantities. In order not to detract the reader from the 
physics, we work with Newtonian definitions here, and the gener¬ 
alization to the relativistic case is discussed in Appendix|A|instead. 
Unless explicitly stated otherwise, the analysis is based on models 
sll.2_2Da and sll.2_3D, i.e. we always refer to model sll.2_2Da 
and not to model si 1.2_2Db when talking about the 2D case. 

5.1 Mass and Enthalpy Flux into the Ejecta Region 

The first quantities to consider are the mass fluxes M in and M out in 
the downflows and outflows, 

Vin/out = I pv r r 2 dn, (5) 

Jv r ^0 

where p is the density, and the less and greater signs refer to down¬ 
flows and outflows respectively. Furthermore, we define total en- 
thalp^and energy fluxes F h and F e , 

Fh,m/om = f [p(e + v 2 /2 + ®) + P] v r r 2 dfi (6) 

Jv r fo 

4 The quantity h to t = e + P/p + v 2 /2 + ®, which we shall usually designate 
in this paper as “total enthalpy” is also referred to as Bernoulli integral or 
stagnation enthalpy in other contexts (including the case without gravity). 
In this paper, we prefer the term “total enthalpy” to keep the terminology 
compact and stress its close relation to the total energy per unit mass. 
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,in/out = [ p(e + v 2 /2 + ®)v r r 2 dn 
0 


(7) 


where 6 is the mass-specific internal energy, v is vectorial fluid ve¬ 
locity, and O is the Newtonian gravitational potential. The rationale 
for including the gravitational potential in these fluxes is that there 
is a conservation law for the total energy density e tot = e + v 2 /2 + ®, 


d_ 

dt 


P\ e 


i + ° 


V- 


P 6 


y +0|V + PV 


= 0 , ( 8 ) 


if the gravitational potential is time-independent and the conversion 
of rest mass energy is either disregarded or nuclear rest masses are 
included in the internal energy. Since the diagnostic explosion en¬ 
ergy is defined as an integral over the total energy of the ejecta, the 
energy budget of the ejecta naturally involves the total enthalpy flux 
from lower regions of the gain layer to the “ejecta region” with e tot 
if the boundary of the ejecta region remains at a roughly constant 
radius. 

In Figure 


12 


we show M out and the average total enthalpy 
and energy h tot and e out (defined as h out = F Kout /M out and e out = 
F e ,out/M ou u and excluding rest mass contributions) in the outflows 
at a radius of 400 km as a function of time. This radius has been 
chosen because recombination into cr-particles, which roughly sets 
the final mass-specific total energy in the ejecta |Scheck et al.| 
|2006j, is already complete at this point, so that F h , out roughly 
represents the rate at which net total energy is pumped into the 
ejecta region assuming steady-state conditions (i.e. small varia¬ 
tions of F hpn with time and radius). This is indeed a very good 
approximation as the comparison for F hout with the time derivative 
dF exp i/d£ of the explosion energy as shown by the middle panel 
in Figure [12] dE cxp \/dt correlates extremely well with F hmi , but 
is slightly smaller. The difference is due to the accumulation of 
shocked material with slightly negative total energy and energy ex¬ 
change with the downflows due to turbulent diffusion. The simi¬ 
larity of dE exp \/dt and F hjOUt also indicates that explosive nuclear 
burning does not play a major role for the 11.2M 0 model in agree¬ 
ment with earlier 2D simulations with the Vertex-CoCoNuT code 
( [Muller, Janka & Marek|20T2] >. 

On average, the total enthalpy flux into the ejecta region is 
larger in 3D than in 2D as expected from the different evolution 
of the explosion energy. Interestingly, the relative difference be¬ 
tween 2D and 3D in the mass outflow rate M out is even larger 
than for F h , ou t- The smaller outflow rate in 2D is partially compen¬ 
sated by a larger average mass-specific total enthalpy and energy 
in the outflows, which can be larger than the recombination energy 
(7... 8.8 MeV/nucleon). Thus, care must be exercised in explain¬ 
ing differences between 2D and 3D based on the mass outflow rate 
or the total mass in the gain region alone (Scheck et al. 2006, [Meh| 
| son, J anka & Marek 2015) assuming the same contribution to the 
explosion energy per unit mass from nucleon recombination into a- 
particles and heavy nuclei irrespective of the dimensionality. That 
assumption would require similar average enthalpies and energies 
in the outflows in 2D and 3D, which is clearly not the case in gen¬ 
eral; the differences can be as large as several MeV/nucleon. Re¬ 
combination still sets the scale for the asymptotic total energy per 
unit mass of neutrino-heated ejecta, but hydrodynamic effects mod¬ 
ify its precise value in 2D and 3D in different directions as we shall 
explain below. 

These differences are all the more astonishing because the 
volume-integrated neutrino heating rate Ghe at in the gain region 
(Figure [13]) is very similar in 2D and 3D. Especially at late times, 
Cheat is consistently higher in 2D than in 3D (as is the time- 
integrated neutrino energy deposition). Assuming that the outflow 


rate is determined by the total heating rate and the binding energy 
£ g ain at the gain radius as M out ~ Cheat /kgainl, the lower outflow rate 
in 2D suggests that the material at the gain radius is more strongly 
bound in this case. This is borne out by Figure [14] which shows 
that the binding energy at the gain radius is larger in 2D by up to a 
factor of ~2 at late times. This is partly due a stronger recession of 
the gain radius r gain (bottom panel of Figure[l4]) as a result of which 
the typical energy scale GM/r gain (M being the neutron star mass) 
at the gain radius is larger. However, it is evident that this effect can¬ 
not fully account for the difference in e gain . The small value of e gain 
in 3D, indicates that it is not determined by the gravitational energy 
scale alone. Indeed, a much better estimate for the scale of e gain can 
be obtained if we suppose that the Bernoulli integral (including rest 
masses) at the gain radius is roughly zero, 

v 2 P 

€ + — + — +0 — 0. (9) 

2 p 

If we split the internal energy e into a thermal component e t herm 
and a rest mass contribution e m (normalized to 56 Fe), and assume 
vanishing velocities as well as an ideal gas equation P - e t herm/3 
of state for photons and non-degenerate relativistic electrons and 
positrons, this leads to, 

P GM 

^therm + ^rm + — — 0? (10) 

P ^gain 

4 GM 

g ^therm + 6rm — _ 0, (H) 

■3 T'gain 


^itherm ■ 





( 12 ) 


where M is the neutron star mass. For an electron fraction of 
Y e = 0.5, 6^ would be identical to the recombination energy of 
e rec ~ 8.8 MeV/nucleon from protons and neutrons with equal mass 
fractions into iron group elements, and for our purposes this still 
provides a sufficient approximation even though we have Y e < 0.5 
at the gain radius. For the binding energy e gain (in which we ex¬ 
cluded rest masses), we thus obtain 

GM ~ 3 GM 

^gain ~ ^itherm — ~ 7 ^rec • D -V 

Y gain 4 4r gain 

As shown in Figure p~4| this still overestimates k gain | a bit, but ac¬ 
counts for the slow rise of k gain | compared to the gravitational en¬ 
ergy scale GM/r gain during the contraction of the neutron star. 

While the different absolute value of the binding energy at 
the gain radius is part of the explanation for the different mass 
outflow rates, there is also an additional effect at play. In gen¬ 
eral, the mass outflow rate will only be approximately given by 
M out ~ Gheat/kgainl, and one can introduce an efficiency parameter 
rj ou t to compare the actual mass outflow rate with this fiducial rate, 


hout 


kgainl^out 

Gheat 


(14) 


h ou t is plotted in Figure[l5] On average, the outflow efficiency rj out is 
also considerably larger in 3D (where it fluctuates around r] out ~ 1) 
than in 2D (r] out & 0.5). 

Our analysis of the outflows has thus revealed two reasons 
for lower explosion energies in 2D: The mass loss rate, (and hence 
the energy flux into the ejecta region) is lower because the ejected 
material is initially bound more tightly at the gain radius before 
being lifted out of the gravitational potential well, and for a given 
binding energy e gain at the gain radius, the conversion of neutrino 
heating into an outflow is less efficient. 
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Figure 16. Snapshots of the specific entropy 5 (left column, measured in ^/nucleon) and the radial velocity v r (right column, measured in cm s -1 ) for model 
si 1.2_2Da (left halves of the individual panels) and for a slice of model si 1.2_3D (right halves) at a post-bounce time of 400 ms. The same data is shown in all 
plots, only the zoom level is different. Note the broad equatorial accretion downflow and the formation of a secondary accretion shock at a radius of -100 km 
in 2D. 


5.2 Causes for Weak Explosions in 2D 

These two effects, as well as the higher asymptotic energy per unit 
mass in 2D can be traced to the constrained axisymmetric flow ge¬ 
ometry and a fundamentally different behaviour of the accretion 
downflows in 2D compared to 3D. The different flow morphology 
is illustrated qualitatively in Figure[l6] which shows snapshots (for 


different spatial scales) of the radial velocity and entropy in 2D and 
3D for a representative post-bounce time of 400 ms. 
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5.2.1 Morphology and Dynamics of Outflows and Accretion 
Downflows in 2D and 3D 

These snapshots reveal that the interface between the outflows 
and the colder, low-entropy becomes turbulent due to the Kelvin- 
Helmholtz instability in 3D, which distorts the downflows as they 
approach the proto-neutron star (as already mentioned in Sec¬ 
tion [3j, whereas Kelvin-Helmholtz instabilities are noticeably ab¬ 
sent at the shear interfaces between the broad equatorial downflow 
and the polar bubbles in 2D (middle and bottom row in Figure |T6j. 
The tendency of the downflows to become more turbulent in 3D 
has been recognized already by |Melson, Janka & Marek| ( [2015[ ) in 
their 9.6M 0 model, although the morphological difference between 
2D and 3D is much more pronounced in a continuously accreting 
model like the 11.2M 0 progenitor. Moreover, although there may 
be a deeper connection between the two phenomena, the stability 
of the shear interfaces in 2D is likely due to a different reason than 
the emergence of large-scale structures in the pre-explosion phase 
( |Hanke et al.|2012j ) that has been explained by the inverse turbulent 
energy cascade in 2D (Kraichnan 1967 ). Despite the inverse turbu¬ 
lent cascade, subsonic shear layers/interfaces remain prone to the 
Kelvin-Helmholtz instability in 2D; and 2D supernova simulations 
are easily able to resolve the instability ( [Muller, Janka & Heger| 
|2012||Fernandez|2015| ) even without extraordinary high resolution. 

The situation changes, however, for the supersonic shear in¬ 
terfaces between the downflows and the neutrino-heated bubbles 
that we encounter during the explosion phase. Here, the classi¬ 
cal growth rate a> = kAu/2 (where k is the wave vector, and A u 
is the transverse velocity jump across the interface) in the vortex 
sheet approximation for incompressible flow is no longer applica¬ 
ble. Instead modes with a sufficiently high effective Mach number 
Ma = Aucos6/c s (where c s is the sound speed and 6 is the angle 
between the wave vector and the vectorial velocity jump) are stabi¬ 
lized ( |Gerwin|1968] ), although the stability analysis is more com¬ 
plicated if finite-width shear layers are considered ( 


I Blumen, Drazin & Billings|1975[|b razin & Dave y|1977[ |Ch oud- 

|hury & Lovelace|l984[|Balsa & Goldstein 1990) ^This implies that 
the Kelvin-Helmholtz instability can be partially or completely sup¬ 
pressed in 2D (where cos 6 =1), while there are always unstable 
modes in 3D since cos 6 can be arbitrarily smallj^Jln principle, it is 
conceivable that numerical diffusivity and viscosity further help to 
suppress the Kelvin-Helmholtz instability more strongly and ear¬ 
lier than the physics might dictate, but the fact that the instabil¬ 
ity evidently operates in the 3D model provides evidence that the 
numerical resolution cannot be faulted for the behavior of the 2D 


5 It is noteworthy that laser-driven plasma experiments may be able to 
capture this effect and quantify the reduced growth or suppression of the 
Kelvin-Helmholtz instability in 2D ([Malamud et al.|2013|. 

6 Loosely speaking, a small value of cos 0 guarantees that sound waves 
on either side of the vortex sheet can propagate in both directions along 
the wave vector k of a given perturbation mode to mediate the pressure 
feedback required for the growth of the Kelvin-Helmholtz instability: For 
a given vectorial velocity ±u/2 of the fluid on either side of the interface, 
the sound waves with direction n and velocity c s in either of the fluids will 
have a velocity component (±u/2 + nc 5 ) • k/|k| = ±w/2cos# + c s n • k/|k| 
along the direction of k in the rest frame. If cos 6 is sufficiently small, this 
velocity component can take on either sign depending on n in both fluids. In 
2D, we always have ±u/2 • k/|k| = u/2, and sound waves cannot propagate 
in both directions any longer for sufficiently large u. Note that this is only a 
heuristic explanation that cannot predict the critical Mach number correctly; 
Section B in |Gerwin]jl968) and the other aforementioned references should 
be consulted for a rigorous derivation of the dispersion relation. 


models. Furthermore, the stability of the accretion downflows is a 
persistent feature even in high-resolution 2D models with continu¬ 
ing accretion ( [Bruenn et al.|2014| >; it is thus without doubt physical 
in origin. 

The “turbulent braking” of the downflows in 3D is reflected 
quantitatively in radial profiles of the average velocity v in / 0 ut, en¬ 
tropy Si n / out and mass-specific total energy Got,rm,in/out of the down¬ 
flows and outflows. We define these quantities as density-weighted 
averages (denoted by bars) of the radial velocity v r , the specific 
entropy s , and the total energy e tot as 


Hn/out — 


^in/out 


Got,in/out — 


X,go^ l ’-' dn 

l $ ,r dn 

l S Qp sdn 
: l, $ oP dQ 

X,§0 P e tovm ^ 
X, ;g0 P 


(15) 


(16) 


(17) 


and show the results for models si 1.2_2Da and si 1.2_3D at a post¬ 
bounce time of 400 ms in Figure^] Moreover, we consider radial 
profiles of the mass and total enthalpy fluxes M in / out and Fh,rm,in/out 
in both streams in Figure [l8j these are computed according to 
Equation {5} and D- However, for computing radial profiles we 
include rest mass contributions in the total energy and the total en¬ 
thalpy flux (as denoted by the additional subscript rm). This defini¬ 
tion has the advantage that both M in>out and F h ,rm,in/out are constant 
in the limit of stationary streams without mass, energy, and mo¬ 
mentum exchange, so that changes in these fluxes serve as useful 
indicators for lateral mixing between the outflows and downflows. 

Due to turbulent braking (i.e. by an effective turbulent eddy 
viscosity), the average infall velocity in the downflows reaches only 
1.4 x 10 9 cm s -1 in 3D, and decreases in magnitude once the down¬ 
flows penetrate further down than a radius of ~200 km. By contrast, 
the downflows reach a sizable fraction of the free-fall velocity in 2D 
before they are abruptly decelerated at a secondary accretion shock 
at r « 100 km. However, the outflow velocities are also higher in 
2D. 

During the phase considered here, the entropy of the down¬ 
flows does not vary considerably in 2D between r ~ 100 km and 
r « 1000 km (where the equatorial downflow forms from two 
converging lateral flows). Similarly, the “flux-averaged enthalpy” 
Fh, rm ,in/out/Mn/out (bottom panel of Figure [?8j does not change ap¬ 
preciably in this region. This is a further indication for a lack of dis¬ 
sipation by turbulent eddy viscosity and of lateral mixing between 
the downflows and outflows. By contrast, the slope in s in and s out 
and the “flux-averaged” total enthalpy F h;rm; j n / out /M in / out points to 
lateral mixing between the downflows and outflows in 3D. The in¬ 
crease of F h?rm ,in/out/M n /out in the downflows during the infall from 
the shock is mirrored by a decrease of F h>m yn/out/Mn/out with radius 
in the ouflows due to turbulent mixing of the hot ejecta with colder 
matter from the downflows: This explains why the neutrino-heated 
ejecta contribute only ~5 ... 6 MeV/nucleon to the diagnostic ex¬ 
plosion energy, instead of the 7... 8.8 MeV/nucleon available from 
nucleon recombination. 

There is thus ample evidence that turbulent viscosity and dif¬ 
fusion brake the accretion funnels and transfer energy from the out¬ 
flows. It is tempting to invoke this as an explanation for the lower 
binding energy at the gain radius in 3D (Figure[l4|. At first glance, 
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the fact that both the downflows and outflows are less strongly 
bound in 3D at the bottom of the gain layer (bottom panel of Fig¬ 
ure EU may seem to conflict with this assumptions, but the lower 
binding energies of the outflows at small radii is to be expected be¬ 
cause the supply for outflow comes from freshly accreted matter 
that has undergone turbulent braking in the downflows. Turbulent 
braking and turbulent diffusion are perfectly acceptable explana¬ 
tions for internal energy distribution within the gain region (but 
not the higher enthalpy flux into the ejecta region, see below), and 
may thus account for the lower |e gain | in 3D and hence for the higher 
mass outflow rates. 

Moreover, the turbulent braking in 3D may have implica¬ 
tions for the final neutron star masses. In Section [373] we assumed 
that the “mass cut” occurs roughly when the shock accelerates the 
newly swept-up material to the escape velocity. Without efficient 
lateral momentum transfer and without pressure support from an 
expanding hot bubble from below, it seems inevitable that material 
over existing downflows must eventually fall onto the neutron star 
if this condition is not met. Since the free-fall time scale at radii 
of several thousands of kilometers (where the initial mass cut esti¬ 
mated in Section |33] is located) is of the order of seconds, accretion 
must necessarily last over a correspondingly long duration. On the 
other hand, if the downflows are braked by a turbulent eddy vis¬ 
cosity below a certain radius in 3D, slowly moving matter in the 
downflows may instead be entrained by the high-entropy bubbles 
in regions where the angular-averaged velocity is already positive, 
so that the residual accretion onto the neutron star may cease ear¬ 
lier and the total mass accreted during the explosion phase may be 
considerably lower than estimated in Section [33] Longer 3D sim¬ 
ulations will be necessary to investigate this hypothesis. 

However, the internal redistribution of energy within the gain 
region in 3D cannot account for the significantly higher total en¬ 
thalpy flux in the outflows and the faster rise of the explosion en¬ 
ergy: The higher outflow rate will come at the expense of energy 
loss from the outflows to the downflows - the overall conserva¬ 
tion law cannot be cheated. Consequently, there must be additional 
mechanisms that remove energy from the gain region in 2D and re¬ 
duce the outflow efficiency r] out (Figure [T5J compared to 3D. The 
different dynamics of the outflows and downflows nonetheless re¬ 
mains a crucial element of the explanation because it provides the 
basis for three mechanisms discussed in the subsequent sections. 

5.2.2 Energy Loss by Wave Excitation 

The lack of turbulent braking in 2D implies that the accretion fun¬ 
nels either hit the neutron star directly with a high impact velocity 
or are decelerated abruptly in a secondary accretion shock (top row 
of Figure [16} . Even in the latter case, thin accretion funnels still 
penetrate the hot, neutrino-heated matter all the way down to the 
gain radius, cutting off a confined high-entropy bubble from the 
outflows. In the snapshots shown in Figures |T6] and [19] (with an 
even higher zoom level), these narrow downflows strike the proto¬ 
neutron star surface with velocities of up to 6 X 10 9 cm s -1 . As a 
result they overshoot considerably into the convectively stable cool¬ 
ing layer and excite strong wave activity. The emission of strong 
acoustic waves that steepen into shocks and then dissipate is imme¬ 
diately evident from the top right panel of Figure |T6] but the decel¬ 
eration of the downflows also excites strong g-modes in the neu¬ 
tron star surface region (a phenomenon that has been thoroughly 
analyzed in the context of gravitational wave emissi on, cp. [M arek, 

Janka & Muller j2009]|Murphy, Ott & Burrows |2009[|Mullei~J anka 

& Marek 2013). In 3D, overshooting is much less pronounced, and 



r [km] 


Figure 17. Radial profiles of the average velocity (top panel), entropy (mid¬ 
dle panel), and total energy (bottom) per nucleon in the outflows (thin lines) 
and downflows (thick lines) in 2D (black) and 3D (red) at a post-bounce 
time of 400 ms. Note that rest mass contributions are included in the total 
energy here. 

so is the excitation of acoustic waves and g-modes. This can be seen 
from the conspicuous absence of sawtooth-like features in the ve¬ 
locity and by considering the radial velocity dispersion ((v r -(v r )) 2 ), 
which is significantly smaller in 3D below the gain radius (top 
panel of Figure [20} . This result is in agreement with other 3D 
simulations of supernova explosions using self-consistent neutrino 
transport (Mels on, Janka & Ma rek 2015 } and pa rameterized neu- 
trino heating ([Murphy, Dolence & Burrows|2013||Handy, Plewa &] 
|Odrzy wolek 2014} and linear theory, which suggests that the exci¬ 

tation of waves (g-modes in particular) at convective boundaries is 
strongly sensitive to the convective Mach number Ma and becomes 
very efficient for Ma ~ 1 ( [Goldreich & Kumar|1990[ |Lecoanet~&| 
|Quataert|2013| >. 

What has been missed so far, however, is that the excitation of 
g-modes constitutes a non-advective energy drain in 2D; it trans¬ 
ports energy from the lower layers of the gain regions deep into the 
cooling region without the need to transport mass. If the g-mode 
energy flux is sufficiently high, it provides a very likely explana¬ 
tion for the permanently higher binding energy at the gain radius in 
2D. Unfortunately, the g-mode energy flux in our simulations can¬ 
not readily be quantified; this would not only require performing 
a full spherical Reynolds decomposition, but also detailed knowl¬ 
edge about the dispersion relation of the g-modes, which is beyond 
the scope of this paper. However, since the transfer of the kinetic 
energy from the downflows into g-modes involves P dV-work onto 
the neutron star surface and any turbulent energy flux into deeper 
layers should show up in correlated pressure and velocity fluctua- 
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Figure 18. Overview of the mass and energy fluxes in the outflows and 
downflows in 2D. The top panel shows the total enthalpy fluxes / 7 A, rm ,i n / 0Ut 
in the outflows (positive values) and downflows (negative values) in 2D 
(black) and 3D (red) at a post-bounce time of 400 ms. The middle panel 
shows the mass inflow/outflow rates Mi n / out , and the bottom panel shows 
the average flux-weighted total enthalpies Fh,rm,in/out /Note that rest 
mass contributions are included in the total entropy here, unlike in Fig¬ 
ures [12] and [T7] 

tions in a transition layer between the convective gain region and 
the convectively stabilized cooling layer, we can formulate a crude 
estimate for the g-mode flux by computing what is nominally an 
acoustic luminosity, namely, 

L Piv = r 2 JsP6v r dQ, (18) 

where 8P and 8v r denote the deviations of the pressure and radial 
velocity from their respective angular averages. The resulting esti¬ 
mates for the flux are shown in the bottom panel of Figure |20| and 
point to a sizable energy flux of the order of several 10 50 erg s -1 
from the vicinity of the gain radius into the deeper layers of the 
proto-neutron star surface (carried by g-modes) and to the outer 
regions of the gain layer (carried by acoustic waves). Such large 
fluxes are comparable to the typical total enthalpy flux in the out¬ 
flows and even to the volume-integrated neutrino heating rate, and 
therefore need to be accounted for in the total energy budget of the 
gain region. 

Incidentally, the excitation of acoustic waves also provides 
an explanation for the high entropy (Figure [IT} and total enthalpy 
(Figures [T2| and p~8]> in the outflows in 2D, which can still increase 
somewhat at radii where neutrino heating is essentially irrelevant. 
The dissipation of the the acoustic waves in the expanding hot bub¬ 
ble helps to increase the total energy and entropy in the ejecta re¬ 
gion beyond the ~7... 8.8 MeV available from nucleon recombi¬ 


nation, but since these waves carry only part of the energy lost due 
to by the downflows due to their interaction with the convective 
boundary this effect cannot compensate for the lower mass outflow 
rate in 2D, and the net effect of wave excitation in 2D remains a 
detrimental one. 

The importance of wave excitation at the convective bound¬ 
ary will inevitably vary between different 2D models depending on 
the explosion geometry and the duration of accretion. If neutrino 
heating is strong, the explosion energy rises steeply, and shock ex¬ 
pansion is fast as in the models of |Bruenn et al. (2013 2014]), it 
may play a less prominent role. The formation of secondary accre¬ 
tion shocks as in our 2D models is likely to increase the energy loss 
by wave excitation tremendously because the efficiency of this pro¬ 
cess also depends on the frequency overlap between the convective 
forcing and the excited modes ( [Goldreich & Kumar|1990[|Lecoanet| 
|& Qu ataert 2013). The formation of a secondary accretion shock 
provides for fluctuations with typical frequencies inversely propor¬ 
tional to the short sound-crossing time-scale (of the order of mil¬ 
liseconds) in the confined bubble. Furthermore, the stochastic forc¬ 
ing of g-modes in 2D by one or two strong downflows is presum¬ 
ably also more efficient than in 3D, where there are more smaller 
and uncorrelated downflows. 

Finally, we comment on similarities and differences between 
g-mode and acoustic wave excitation between our models and 
the acoustically-driven explosion models of |Burrows et al.|f2006[ 
2007), where a strong flux acoustic waves, excited by an l - 1 
core g-mode with amplitudes of several kilometers, is responsible 
for shock expansion in the first place. Our 2D models are similar to 
those of Burrows et al. ( 2006 2007| only in the sense that energy 
deposition by acoustic waves contributes to the growth of the explo¬ 
sion energy, but different from their simulations the acoustic con¬ 
tributions remains subdominant compared to the volume-integrated 
neutrino heating rate, which is more than four times larger at the 
time shown in Figure |20] (a constellation that Burrows et al. |2007| 
anticipated when postulating a “hybrid mechanism” with combined 
heating by neutrinos and acoustic waves). Moreover, the excitation 
mechanism for acoustic waves is genuinely different in our case; 
they are excited directly by the interaction of the downflows with 
the convectively stable neutron star surface layer without the need 
to channel the accretion power through an £ = 1 core g-mode as a 
“transducer” as in the models of Burrows et al. ( 2006 2007). Such 
a large-amplitude core g-mode is not found in our simulations (and 
could not have arisen simply because of the spherically symmet¬ 
ric treatment of the neutron star core), and the outer g-modes of 
rather modest amplitude excited in our models could not act as an 
efficient transducer in the vein of |Burrows et a l. (2006 2007) due 
to neutrino losses (see below). It is interesting to note that direct 
excitation at the convective boundary still allows acoustic waves to 
contribute (albeit at a minor level) to the explosion energy in 2D 
even without such a transducer. 

Acoustic energy deposition also remains a secondary effect 
insofar as this direct excitation mechanism works efficiently only 
after the onset of the explosion once the typical Mach number of 
the downflows is sufficiently high. Moreover, contrary to Burrows 
|et al. 1 (2006 2007) the net effect of wave excitation in our models is 
still harmful because the power pumped into outer g-modes consti¬ 
tutes an energy drain that outweighs the rate of energy deposition 
by acoustic waves by far. Different from their models where the 
energy in the core g-mode is eventually “recycled” into an acoustic 
energy flux that drives shock expansion, the energy pumped into the 
outer g-mode is manifestly lost due to neutrino cooling in our case, 
and the mode coupling analysis of | Weinberg & Quataert| ( [2008] ) 
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suggests that this should also happen if the core g-mode excited 
due to non-linear mode coupling. 


5.2.3 Steric Hindrances 

In addition to energy loss by wave excitation, which contributes to 
the higher binding energy at the gain radius, the growth of the ex¬ 
plosion energy in 2D is further hampered by the fact that much of 
the neutrino energy deposition occurs in regions where the heated 
matter cannot directly escape in an outflow, i.e. either directly in 
the accretion funnels or in high-entropy bubbles confined by down¬ 
flows and a secondary accretion shock like the equatorial bubble in 
Figures [Tb] and [T9| a phenomenon for which we borrow the term 
“steric hindrance” from chemistry. In principle, such bubbles could 
eventually push the secondary accretion shock out by undergoing 
“secondary shock revival”, but as long as the amount of heating is 
insufficient, the bubble cannot break through the surrounding and 
overlying downflows. 

In the snapshot shown in the right panel of [19] the surface 
fraction covered by the confined bubble and the downflows exceeds 
50%, and the heating rate per unit mass is also largest in the down¬ 
flows. Since the fast downflows generally occupy a surface fraction 
of >50% in 2D outside the typical location of a secondary shock 
at < 100km (see Figure [2T]>, roughly half of the neutrino heating 
is not used to power outflows in 2D, and consequently the out¬ 
flow efficiency oscillates around 77 ou t ~ 0.5 with some excursions 
to higher values during the early explosion phase (Figure [HJ. By 
contrast, the neutrino-heated material can escape unhindered in any 
direction in 3D apart from some limited turbulent energy and mo¬ 
mentum loss to the downflows, and the resulting outflow efficiency 
is of order 77 out ~ 1. 

5.2.4 Constriction of Outflows and Vertical Mixing 

Finally, the outflows in axisymmetric 2D simulations are less “sta¬ 
ble” than in 3D in other respects as illustrated in Figure [22] While 
the Kelvin-Helmholtz instability between the downflows and out¬ 
flows is largely suppressed in 2D, this also implies that plumes 
of cold material can penetrate far into the neutrino-heated high- 
entropy bubbles provided that they develop in the first place. Be¬ 
cause of the symmetry of the system, these plumes are actually 
toroidal structures, and can therefore completely constrict an out¬ 
flow if they reach the symmetry axis (Figure|22]>. Similarly, a down¬ 
flow that wanders toward the pole can also constrict a polar outflow 
and cut it off from fresh supply of neutrino-heated material. 

These events typically reduce the surface fraction covered by 
outflows at the recombination radius (where they start to contribute 
to the diagnostics explosion energy) for a considerable amount of 
time and thus reduce the rate of increase of £ exp i. Very often the 
explosion geometry changes dramatically after such an event and 
the surface fraction of the outflows remains small permanently. In 
some cases, a high-entropy bubbles is cut off completely from the 
supply of neutrino-heated matter from below (Figure [23]) for sev¬ 
eral seconds. Even if an outflow is eventually reestablished in the 
same direction, or if it is strong enough to survive because the cold 
plumes reach the axis at a relatively large radius (Figure [22]), the 
ejecta will then typically contain a large amount of cold material 
whose total net energy is barely positive, and the growth of the ex¬ 
plosion energy will still be delayed. Such events explain excursions 
or even a permanent drop of the the average total enthalpy h tot in 
the outflows to low values < 0.3 in 2D (bottom panel of Figure \V2). 


In 3D, the lack of symmetry as well as the Kelvin-Helmholtz 
instability prevent the constriction of outflows by cold plumes. 
While the Kelvin-Helmholtz instability provides for some level of 
energy and momentum exchange between the accretion funnels and 
the expanding high-entropy bubbles as discussed in Section |5~.2.1[ 
it also prevents cold plumes from penetrating overly far into the 
neutrino-heated bubbles. 

Mixing between downflows and outflows is thus not com¬ 
pletely absent in 2D, it merely takes on a different guise and occurs 
only episodically, but with a more catastrophic effect than in 3D. 
Interestingly, there even appears to be an effect that compensates 
somewhat for the suppression of the Kelvin-Helmholtz instabilities 
in 2D due to the supersonic velocities in the downflows: Rayleigh- 
Taylor instabilities between the high-entropy bubbles and the cold 
overlying post-shock matter develop more readily in 2D. This is a 
natural consequence of higher entropies in the neutrino-heated bub¬ 
bles in 2D (middle panel of Figure [17]), which implies a higher At¬ 
wood number between the bubbles and the colder post-shock mat¬ 
ter. Thus, the lack of mixing by Kelvin-Helmholtz instabilities in 
2D and the entropy boost due the dissipation of acoustic waves can 
also have a detrimental side effect on the robustness of the explo¬ 
sion. 

5.2.5 Absence of a Spherically Symmetric Neutrino-Driven Wind 
in 2D 

In the most extreme cases of outflow constriction in 2D, the out¬ 
flows are shut off altogether, and the outflow surface fraction drops 
to zero permanently, or at least over several seconds (bottom panel 
of Figure [2l] and Figure [23]). This does not imply, however, that 
the explosion has failed; it only implies that neutrino heating is 
not strong enough to establish a wind that prevents the fallback of 
slowly moving matter in the wake of the shock. The pockets of 
cold, slowly moving matter from the C/O shell in the 2D models 
that will undergo this kind of “early fallback” (bottom right panel 
of Figure [23]) only contain a few hundredths of a solar mass by the 
end of the simulations, and therefore will not change the neutron 
star mass considerably. Moreover, the mass accretion rate onto the 
secondary accretion shock is so low at late times that it can start 
to expand after “secondary shock revival”, thus re-establishing an 
outflow (bottom right panel of Figure [23]). 

While not indicative of a failure of the explosion, the small 
or vanishing outflow surface fraction in the long-time simulations 
nonetheless indicates (like the models of [Bruenn et al.|20l4] ) that 
the separation of outgoing and infalling mass shells in 2D works 
differently from the usual picture where a high-entropy neutrino- 
driven wind with an approximately spherical flow geometry even¬ 
tually develops. The polar outflows can be viewed as a confined 
wind driven jointly by neutrino heating and acoustic waves, but 
they never cover the entire sphere, and because of their flow geom¬ 
etry and the strong activity of acoustic waves, the outflow dynamics 
is completely different from spherical winds driven purely by neu¬ 
trino heating. 

From the foregoing, it is clear that the inhibition of the 
neutrino-driven wind in 2D is probably largely artificial; and we 
only mention this peculiarity for that reason. As discussed in Sec¬ 
tion [5T27T] the presence of a larger effective eddy viscosity in 3D 
could terminate accretion earlier than in 2D (where supersonically 
infalling matter is hardly decelerated by lateral momentum trans¬ 
fer), or at least decelarate infalling matter sufficiently to be swept 
along by an incipient spherical wind after a few seconds. Moreover, 
our models likely underestimate the diffusive neutrino luminosity 
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Figure 20. Top: Radial velocity dispersion <(v r - (v r » 2 ) in 2D and 3D at 
a post-bounce time of 400 ms. Bottom: Radial profiles of the “acoustic” 
energy flux r 2 f 6P6v r dD in 2D and 3D at a post-bounce time of 400 ms. 
The curves show temporal averages over several time steps. 

from the neutron star core and hence the neutrino heating at late 
times because we ignore the effect of nucleon correlations |Bur-| 
|rows & Sawyer|1998||1999|[Reddy et al.|1999| ), which shorten the 
proto-neutron star cooling time-scale considerably ( Hudepohl et al.| 
12009] ). It is conceivable that the concomitant increase of the wind 
mass loss rate could still lead to a volume-filling outflow for more 
realistic neutrino opacities after a few seconds even in 2D. 

5.2.6 Reduced Cooling due to 3D Turbulence? 

Based on a successful supernova simulation of a 9.6M 0 star, |Mel-| 
| son, Janka & M arek ( 2015) recently suggested that the more ef¬ 
ficient braking of the accretion downflows can be responsible for 
slightly higher explosion energies in 3D because the less violent 
impact of the downflows on the neutron star surface lead to reduced 
cooling. In our comparison of models sll.2_2Da and sll.2_3D we 
observe some of the same symptoms noted by these authors, i.e. tur¬ 
bulent braking of the downflows and a reduced cooling rate gcooi at 
late times (Figure [T3J. Obviously, this raises the question whether 
the mechanism proposed by |Melson, Janka & Marek] ( |2015| ) also 
operates in our 3D simulation, and how the physical processes we 
discussed so far in Sections r5.2.1l to l5'.2.4l are related to it. Unfor¬ 
tunately, a comparison with |Melson, Janka & Marek| ( |2015| ) is not 
straightforward. While they found a sizable increase of the explo¬ 
sion energy of 10% in 3D compared to 2D, their explanation in¬ 
volved relatively tiny differences in some quantities (e.g. the gain 




Figure 21. Top: Surface fraction occupied by the outflows in models 
sll.2_2Da (black) and sll.2_3D (red) at a radius of 400 km. The surface 
fraction is relatively stable with some fluctuations around 0.5 in 3D. In 
2D, it reaches similar values while the outflows are stable, but occasionally 
drops to significantly smaller values as a result of outflow constriction. Bot¬ 
tom: Long-time evolution of the outflow surface fraction for the 2D models 
sll.0_2D, sll.2b_2D, sll.4_2D, and sll.6_2D. 


radius and the temperature profiles in 2D and 3D) that cannot be 
confidently diagnosed in simulations like ours where the 2D and 
3D models start to deviate from each other already shortly after 
bounce once prompt convection develops (which was not simu¬ 
lated by |Melson, Janka & Marek|2015) . Nonetheless, there is suf¬ 
ficient evidence that we observe some rather different phenomena 
than |Melson, Janka & Marek| ( |2015| ). 

Essentially, the mechanism proposed by Melson, Janka & 
|M arek| < |2015| ) involves a recession of the gain radius in 3D com¬ 
pared to 2D due to reduced cooling to eject slightly more mate¬ 
rial in the explosion. Our simulations agree with|Melson, Janka &| 
[M arek| ( |2015] ) in showing a smaller volume-integrated cooling rate 
in 3D in the long term as accretion slowly subsides (Figure [T3). 

However, we do no find a faster recession of the gain radius in 
3D during the explosion phase (middle panel of Figure [M}, and the 
situation is ambiguous for the temperature at the gain radius r gain 
(bottom panel of Figure [14]). In 3D, the temperature r gain stagnates 
and falls below the 2D value around 250 ms at a time when the 
explosion is already considerably more vigorous in 3D than in 2D. 
We believe that the stagnation of r gain is more indicative of the 
slower recession of the gain radius rather than of a higher cooling 
efficiency: The higher values of the Bernoulli integral and the total 
energy of the downflows at the gain radius in 3D (Figure [TT] ancfl8]> 


© 0000 RAS, MNRAS 000, 000-000 











































































Dynamics of Neutrino-Driven Supernova Explosions 19 



x [km] x [km] 

Figure 19. Left: Entropy s in kb /nucleon in the vicinity of the proto-neutron star in 2D (left half of panel) and 3D (right half of panel) at a post-bounce time of 
400 ms (identical to Figure Q6] except for the zoom level). Right: Heating/cooling rate in MeV/nucleon in 2D (left half of panel) and 3D (right half of panel). 
Iso-velocity contours for a radial velocity of v r = -10 9 cm s -1 are shown to indicate the location of the accretion downflows. Note that much of the neutrino 
heating occurs in the downflows and the confined high-entropy bubble in the equatorial region and hence does not drive an outflow. 


imply that there is actually more energy per unit mass available 
that can be radiated away in neutrinos as the accreted matter settles 
down in the cooling region. 

Instead, the faster decline of the accretion rate onto the proto¬ 
neutron star in 3D is the dominant factor behind the lower cooling 
rate, making the lower cooling rate a symptom rather than a cause 
of the more vigorous explosion. Detailed comparisons would be 
required to check whether this true for the 9.6M© model of |Mel-| 
| son, Janka & M arek (2015) as well. Since outflow constriction is 
unlikely to happen for a model with robust shock expansion, the 
2D/3D differences found by |Melson, Janka & Marek |(|2015|) a s well 
as in the parameterize simulations of Handy, Plewa & OdrzywoIek| 
( |2014| ) are probably most closely related to the different outflow ef¬ 
ficiency in 2D and 3D, i.e. a more efficient “rerouting” of freshly 
accreted matter into outflows. This tallies with their finding of a 
smaller surface filling factor of the downflows in 3D, which im¬ 
plies that a smaller fraction of the neutrino heating is wasted in 
regions where it cannot directly power an outflow. It also accounts 
for reduced mass accretion into the gain region and hence a reces¬ 
sion of the gain radius in mass coordinate. While this mechanism 
is similar to the one discussed in Section [5'.2.3| for our models, the 
effect is apparently smaller in the simulations of |Melson, Janka &\ 
Marek ( 2015) because the accretion subsides fast enough to avoid 
the formation of secondary shocks and confined high-entropy bub¬ 
bles in 2D, which can reduce the outflow efficiency by a factor of 
~2 in 2D. 


Potentially, wave excitation at the convective boundary could 
also contribute to the 2D/3D differences in the simulations of [Mel- 1 
| son, Janka & M arek (2015). While they take reduced convective 
overshooting in 3D as an indication for reduced wave excitation, 
the effect probably plays a minor role in their case. The rela¬ 
tively small average speeds of the downflows at the gain radius and 
(~10 8 cm s -1 compared to ~10 9 cm s -1 in our model) are bound 
to make the excitation of g-modes rather inefficient and thus rule 
them out as a major energy drain on the gain region in 2D. This is 


also suggested by the fact they find similar internal energies (and 
hence binding energies) in the outer regions of the cooling layer in 
2D and 3D despite the stronger recession of the gain radius in 3D, 
which is quite different from what we discussed in Sections |5.2.1| 
and !5.2.2l 


6 SUMMARY AND CONCLUSIONS 


We have presented a successful 3D GR simulation of the explosion 
of an 11.2M© star using the FMT multi-group transport scheme 
of |Muller & Janka| ( |2015| >. The model has been evolved to almost 
1 s after bounce, and has reached a diagnostic explosion energy of 
1.3 X 10 5 ° erg at that point, which is still growing by the end of the 
simulation. The baryonic neutron star mass M by at the end of the 
simulations has reached 1.33M© and estimates of the final neutron 
star mass yield M by « 1.41... 1.48M© and a gravitational mass not 
exceeding 1.34M©, which is compatible with the measured neutron 
star mass distribution ( [Schwab, Podsiadlowski & Rappaport|2010| . 
The fact that we obtain an explosion for this progenitor with a rela¬ 
tively accurate multi-group transport scheme further illustrates that 
even the non-exploding state-of-the-art models ( [Hanke et al.|2013[ 
|Tamborra et al.|2014a|b| > with the best available neutrino transport 
and microphysics are apparently very close to shock revival, some¬ 
thing which is also suggested by the recent successful 3D explosion 
models of the Garching (Melson, Janka & Marek 2015, Melson 


|et al.[2015| ) and Oakridge groups ( [Lentz et al.|2015 


A comparison of the explosion dynamics after shock revival 
with 2D long-time simulations of different progenitors with ZAMS 
masses between 11.0M© and 11.6M© revealed a faster and more 
stable growth of the explosion energy in 3D compared to 2D. Be¬ 
cause accretion downflows and neutrino-driven outflows coexist 
over several seconds in 2D, the explosion energy in the 2D mod¬ 
els can still reach values of up to 2 x 10 5 ° erg, but this comes at the 
expense of high neutron star masses (M by > 1.6M©) that are likely 
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Figure 22. Constriction and partial shredding of an outflow in model sll.2_2Da, shown by snapshots of the entropy at post-bounce times of 592 ms, 611 ms, 
628 ms, and 655 ms. A downflow (cyan) originating from a Rayleigh-Taylor plume of cold matter penetrates the hot-neutrino heated bubble in the northern 
hemispheres (top left), constricts the neutrino-heated bubble to a tenuous outflow as it approaches the axis (top right), and eventually a considerable amount of 
cold material is mixed into the outflow (bottom left). While the ejection of matter continues (bottom right), the mixing event lowers the average total energy 
per unit mass in the ejecta. 


incompatible with the observed neutron star mass distribution. A 
detailed comparison of the 2D and 3D models unearthed several 
physical mechanisms responsible for the more robust rise of the 
explosion energy in 3D and the slower growth of the proto-neutron 
star mass, which is a symptom of the faster subsidence of accre¬ 
tion. The specific effects that we find are summarized below, and 
we also provide a schematic visualization of the different flow ge¬ 
ometry and the energy budget between the downflows, the outflows 
and the gain region in Figure [24]to aid the reader’s understanding: 

1. In 2D, the interfaces between the accretion funnels and the 
neutrino-heated bubbles tend to become laminar after shock re¬ 
vival, while they are corrugated by the Kelvin-Helmholtz instability 
in 3D. We ascribe the different behavior to the suppression of the 
purely two-dimensional modes of the Kelvin-Helmholtz instabil¬ 
ity in the supersonic regime ( |Gerwin|1968| ). The effect is thus dis¬ 
tinct from the inverse turbulent energy cascade in 2D (Kraichnan 
|1967| , which has been invoked as an explanation for the different 
behaviour of 2D and 3D models prior to shock revival, since the 


different energy cascade in 2D and 3D is not related to the Mach 
number of the flow. 

2. As a consequence, the effective eddy viscosity and diffusiv- 
ity between the downflows and outflows is larger in 3D than in 2D 
during most phases, i.e. there is more exchange of energy and mo¬ 
mentum between the outflows and downflows. On the one hand, 
this implies that the outflows contribute only ~6 MeV/nucleon 
to the explosion energy in 3D, as some of the net total (i.e ther- 
mal+kinetic+potential) energy gained from nucleon recombination 
of ~8.8 MeV/nucleon is lost to the downflows by turbulent dif¬ 
fusion. On the other hand, the turbulence effectively “brakes” the 
downflows, and they arrive at the gain radius with smaller velocities 
but higher total energy per unit mass than in 2D. 

3. The higher impact velocities of the downflows and the for¬ 
mation of secondary accretion shocks at small radii in 2D lead to a 
more efficient excitation of g-modes and acoustic waves at the gain 
radius that transport energy into deeper regions of the cooling layer 
and into the neutrino-heated ejecta respectively. Our analysis sug- 
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Figure 23. Snapshots of the radial velocity in units of 10 9 cm (left half of panels, lower color bar) and the specific entropy 5 in /^/nucleon (right half of 
panels, upper color bar) depicting the separation of the high-entropy outflow from the gain region in model sll.6_2D and the re-establishment of an outflow 
after “secondary shock revival” in model sll.0_2D. Red isovelocity contours are used to separate outward-moving matter with radial velocities larger than 
10 7 cm s -1 from infalling matter. At 4.4 s (top left), matter is still being ejected in the southern hemisphere, but at 4.5 s (top right) the outflow has become 
extremely thin, and matter in its wake starts falling back onto the neutron star. By 7.5 s (bottom left) a new outflow has developed in the northern hemisphere, 
and the expansion of the secondary accretion shock in the southern hemisphere stops further fallback from the outflow that was cut off earlier. By the end of 
the simulation at 8.2 s (bottom right), the newly formed bubble has expanded further to several thousands of kilometers in diameter. The post-shock velocities 
have become positive in all directions at this point, and only 0.035M o in the downflows are still falling toward the proto-neutron star. 


gests that the energy loss from the gain region by wave excitation 
becomes comparable to the volume-integrated neutrino heating rate 
at late times, and by increasing the absolute value of the binding 
energy \e tot \ at the gain radius significantly reduce the mass outflow 
rate that can be sustained by neutrino heating. In 3D, the turbulent 
energy flux into the gain region is small, and the binding energy 
at the gain radius is smaller by factor of > 2 at late times, which 
allows for a higher mass outflow rate than in 2D. The dissipation 
of acoustic waves in the outflows in 2D provides only for a partial 
“recycling” of the energy lost by wave excitation, but can increase 
the total energy per unit mass in the outflows to values larger than 
the recombination energy of 8.8 MeV/nucleon. 

4. In addition, the outflow efficiency rj 0 ut = M ou t/(Gheat/k g ainl), 


is also higher in 3D (77 out ~ 1 with strong fluctuations) than in 2D 
(?7out ^ 0-5 at late times), i.e. for a given amount of neutrino heat¬ 
ing and a given binding energy at the gain radius, more mass is 
channelled into outflows and contributes to the explosion energy 
in 3D. The low outflow efficiency in 2D stems from the large sur¬ 
face fraction occupied by fast downflows and “confined bubbles” 
between downflows whose expansion is inhibited by the formation 
of secondary accretion shocks. 

5. Episodic mixing between the outflows and downflows still 
occurs in 2D, e.g. by the formation of new downflows as a result 
of the Rayleigh-Taylor between the cold shocked material and the 
neutrino-heated high entropy bubbles. While mixing only occurs 
sporadically in 2D, the consequences of these mixing events are 
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excitation region 

Figure 24. Sketch of the different energy budget between the outflows (yellow), the downflows (red) and the cooling region (blue) in 2D (left) and 3D (right). 
In 3D, turbulent eddy diffusivity leads to a persistent, small energy flux (short solid arrow) from the neutrino-heated outflows to the downflows, whereas 
mixing between the outflows and downflows only occurs episodically in 2D, but has more dramatic consequences in this case because it leads to large-scale 
mixing of cold matter into the outflows (long dotted arrow). There is a considerable energy transfer from the gain region to the cooling region due to wave 
excitation in 2D and an indirect transfer of energy from the downflows to the gain region by the excitation of acoustic waves in 2D (which can lead to larger 
asymptotic energies per unit mass in the ejecta); this is absent in 3D. Moreover, a considerable amount of neutrino heating (red arrows) is wasted in 2D because 
it is deposited in confined bubbles, whereas almost the entire neutrino heating in 3D is used to lift matter out of the gravitational well. 


more catastrophic than in 3D. Not only do they slow down the rise 
of the explosion energy by mixing cold material into the outflows; 
the penetration of accretion funnels into the high-entropy bubbles 
can also lead to the constriction of outflows, sometimes shutting 
them off completely and permanently decreasing the outflow sur¬ 
face fraction to values of < 0.3. 

Our simulations thus provide ample evidence that 3D effects 
can play a beneficial role in core-collapse supernova explosions af¬ 
ter shock revival. However, since our current 3D model has only 
been evolved to ~1 s after bounce and does not yet permit us to de¬ 
duce the final explosion and remnant properties directly because of 
continuing accretion (and forced us to resort to indirect arguments 
about the final neutron star masses), a number of open questions 
remain and invite speculation. Moreover, limited conclusions can 
be drawn from a single 3D simulation of one progenitor. Given the 
recent progress on other fronts in supernova theory, the questions 
and perspectives for future research on the role of 3D effects during 
the explosion phase can be summarized as follows: 

1. Longer 3D simulations with higher resolution are necessary 
to determine final explosion energies, Nickel masses and neutron 
star masses precisely for comparison with observations without re¬ 
course to indirect methods. Our estimate for the final baryonic neu¬ 
tron star mass of 1.48M 0 for the 3D model of the 11.2M 0 pro¬ 
genitor still assumes the accretion of an additional 0.15M o solar 
masses, which is much more than a visual inspection of Figure [To| 
suggests given the very slow rise of the neutron star mass in 3D. If 
the turbulent braking of the downflows terminates accretion before 
the post-shock velocity equals the escape velocity as speculated 
in Section [5.2. 1[ the final baryonic and gravitational neutron star 
mass might be as low as ~1.35M 0 and ~1.24M 0 , respectively. This 
would indicate that a plausible distribution of neutron star masses 
spanning the entire range of observed values down to the lower end 
is within reach of modern multi-D simulations of neutrino-driven 
supernovae. 

2. A more rigorous analysis of the turbulent multi-dimensional 
flow in the spirit of Reynolds decomposition would be highly de¬ 
sirable in order to further bolster our qualitative interpretation of 
3D effects in the post-explosion phase. Such quantitative analysis 
methods have considerably advanced our understanding of the tur¬ 


bulent flow during the accretion phase ( [Murphy & Meakin|2011) . 
After shock revival, the nonstationarity of the flow presents a chal¬ 
lenge for such methods, however. Our relatively crude two-stream 
analysis based on a separation of the outflows and downflows could 
also be improved in order to account more directly and rigorously 
for the turbulent exchange of mass, momentum and energy between 
the two streams, but such an analysis faces a major challenge in the 
form of the complicated flow geometry. 

3. Whether and to what extent the positive 3D effects described 
in this paper come into play obviously depends on whether shock 
revival can be accomplished in 3D in the first place and on the de¬ 
lay compared to the 2D case. If there is a significant delay in shock 
revival, 3D models may not be able to equalize the “head start” of 
the 2D models at least of relatively powerful explosions where the 
diagnostic energy shows first signs of levelling off after ~300 ms 
or less (Bruenn et al. 2014, Pan et al. 2015}. Even in this case, the 
mechanism discussed in this paper could nonetheless help to mit¬ 
igate the “penalty” incurred by the delay of the explosion in 3D 
and allow the models to remain compatible with observational con¬ 
straints. Moreover, if accretion lasts longer — as in the 2D simu¬ 
lations of|Muller, Janka & Marek|(|2012};|Muller, Janka & Heger| 
pOT2t ; [Jmka~ (2012) — the beneficial 3D effects in the phase 
after shock revival may outweigh the “penalty” of delayed shock 
revival. Furthermore, it is conceivable that the problem of missing 
or delayed explosion in 3D may yet be resolved by the inclusion of 
better, multi-dimensional progenitor models with large-scale initial 
perturbations that aid shock revival ( Couch & Ott 2013; Muller & 


Janka|2015l|Couch & Ott|20l5}, unknown microphysics ([Melson 


et al. |2015| ; and strongly SASI-dominated models may even ex¬ 
plode easier in 3D (Fernandez 2 015| ). 

4. The robustness of the mechanisms described in this paper 
needs to be studied further for a wider range of progenitors. The 
2D simulations ( [Burns et al.|2006a[|Muller, Janka & Marek|2012] > 
of the 11.2M 0 progenitor considered here have been particularly 
noteworthy examples for suspiciously low explosion energies and 
long-lasting accretion. This behavior is due to the specific charac¬ 
teristics of progenitors around 11M 0 , including a relatively small 
silicon core and a very pronounced density jump at the Si/SiO in¬ 
terface. These properties result in a small proto-neutron star mass 
M immediately after shock revival and hence low neutrino energies 
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(cp. the scaling of the electron antineutrino energy with M found 
by |Miiller & Janka|2014[ ) as well as a small accretion luminosity. 
Both of these factors contribute to relatively weak neutrino heating 
after shock revival and a small mass outflow rate. The tepid nature 
of our 2D explosions may thus hinge very much on the peculiar 
structure of low-mass supernova progenitors. 

It therefore remains to be seen whether 3D effects provide a sim¬ 
ilarly strong boost for the growth of the explosion energy in other 
progenitors. Whenever 2D models develop the characteristic broad 
downflows and secondary shocks indicative of long-lasting accre¬ 
tion during a relatively weak explosions, such as the 15M 0 and 
27 M q models of [Muller, Janka & Marek| ( |2012| ) and |Janka et al.[ 
( |2012| ), the physical mechanisms identified here likely come into 
play eventually. On the other hand, they may play a negligible role 
if the volume fraction of the downflows drops very quickly as in the 
models of [Bruenn et ah ( 2014| ) an d|Lentz et al.[([2015j) o r the pa¬ 
rameterized models of |Handy, Plewa & O drzywolek (2014 ). Since 
2D supernova simulations of different groups have not yet con¬ 
verged sufficiently to decide whether there is a generic problem 
of “weak explosions” in 2D, it is impossible to judge the generic 
character of our findings. By the same token, however, it cannot 
be ruled out that 2D models should be generically underenergetic 
and overestimate the amount of accretion after shock revival due to 
the mechanisms we identified, and that realistic explosion enegies 
and remnant masses will only be obtained in 3D. If so, prematurely 
confronting the 2D models with the observational constraints could 
lead to wrong conclusions. 


Thus, more work is necessary to substantiate the intriguing per¬ 
spective that 3D effects could help to achieve agreement between 
simulations of neutrino-driven supernovae and observational con¬ 
straints such as explosion energies and neutron star masses. Far 
from offering a complete solution due to the limitations of compu¬ 
tational resources that have always plagued supernova theory, our 
present study can only take a first step in this direction and adum¬ 
brate some of the physical mechanisms that could help to boost 
3D explosions after shock revival. Nonetheless, even our current 
results already serve an antidote against undue pessimism after ini¬ 
tial setbacks in 3D multi-group neutrino hydrodynamics simula¬ 
tions. Along with the recent successful explosions in first-principle 
models, the identification of other beneficial effects of the third di¬ 
mension on the explosion threshold and energetics, and plausible 
ideas for solving the problem of missing explosions with the help 
of multi-D progenitor models and/or non-standard microphysics, 
they are another piece that fits well into the overall puzzle, sug¬ 
gesting that a solution for the supernova problem is slowly taking 
shape. 
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APPENDIX A: DEFINITION OF ENERGIES AND 
ENERGY FLUXES IN GENERAL RELATIVITY 


While the computation of mass fluxes and spherically-averaged 
profiles of hydrodynamic quantities can be readily generalized to 
the relativistic case just by including the correct three-volume ele¬ 
ment, the definition of energies and energy fluxes in the relativistic 
case is less straightforward, and is therefore briefly expounded in 
this Appendix. We use geometric units (G - c - 1) throughout this 
section. 

After adopting a 3 + 1 foliation of space-time and projecting 
the components of the stress-energy tensor into components orthog¬ 
onal and parallel to the 3-hypersurfaces, the energy equation in gen¬ 
eral relativistic hydrodynamics can be written in the formulation of 
|Banyuls et al.| ( |1997] > in terms of a new conserved variable r as 


dyfyr 

dt 


d ( T ^ + Pv 1 ) 


dx l 


= a y/~8\T t ‘° 


d In a 
~dxN 


(Al) 


Here, r is defined in terms of the baryonic mass density p in the 
fluid frame, the Lorentz factor W, the internal energy density e (in¬ 
cluding all rest-mass contributions), the pressure P, and the rela¬ 
tivistic enthalpy h = 1 + e + P/p as 

T = phW 2 -P-pW=p(l+e + P/p)W 2 -P-pW. (A2) 


Furthermore, y and g are the determinants of the three- and four- 
metric, respectively, and the advection term contains v l = v l - /3 1 /a 
instead of the Eulerian three-velocity v l . By pushing the lapse func¬ 
tion a into the temporal and spatial derivatives, it is possible to for¬ 
mulate a strict conservation law (without a source term) analogous 
to in the limit of a stationary space-time with a zero shift vector 
(cp. Equation A35 |Miiller, Janka & Dimmelmeier|2010] where the 
right-hand side reduces to zero in this limit): 

[ sjya (r + D) - sJyD \+[ sf^g [arv 1 + aDv 1 - Dv l + oTV)] = 

(A3) 

Here, we have introduced the baryonic mass density in the Eulerian 
frame D = pW to simplify the equation. 

This suggests that in the limit of a vanishing shift vector and a 
stationary metric, the Newtonian expression e tot = e + v 2 /2 + O for 
the total energy per unit mass (including rest-mass contributions) 
can be generalized to 

(XT 

^tot,rm — + (^ — 1), (A4) 

and the role of the total Newtonian enthalpy in the flux is taken by 


^tot,rm = aT/D + (a - 1) + aP/D. (A5) 

It is noteworthy that the internal energy and rest-mass contributions 
(which enter the equations through r) always appears in conjunc¬ 
tion with factors W and a. Strictly speaking, it is therefore no longer 
possible to formulate the energy equation in general relativity with¬ 
out including rest-mass contributions in the conserved quantities 
and the fluxes by pushing them into a nuclear source term instead 
(as least not in a simple form). Computing fluxes and total ener¬ 
gies excluding rest-mass contributions is therefore somewhat less 
meaningful in general relativity. However, since we have a « 1 and 
v < 0.3c in the gain region, the higher-order relativistic corrections 
are small enough to be neglected, it is still reasonable to compute 
total energies and enthalpies without rest-mass contributions by us¬ 
ing just the thermal energy contribution e t h e rm to the internal energy 
density e instead of e = e t h e rm + e rm . 

For the other quantities, considered in this paper, the gener¬ 
alization is trivial. Mass fluxes through the surface of a sphere or 
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parts of it are computed as 

M= J aDv r r 2 (/> 4 tin, (A6) 

where v r is the radial velocity component measured in the orthonor- 
malised Eulerian frame and 0 is the conformal factor in the xCFC 
metric; and the computation of energy/enthalpy fluxes works anal¬ 
ogously. The density-weighted spherical average X of a quantity^ 
is computed as 


f DX<t> 4 d n 
/Dip 4 dfi ' 


(A7) 
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